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ABSTRACT 


This  thesis  focuses  on  producing  realistic  time-space-position  information 
(TSPI)  data  for  moving  targets  as  inputs  for  a  weapon  flyout  model  to  compute 
the  weapon  Circular  Error  Probable  (CEP)  and  the  probability  of  damage.  The 
velocity  profile  of  the  target  is  modeled  based  on  the  kinematic  constraints  for  the 
type  of  vehicle  and  the  type  of  path  on  which  it  is  traveling.  The  discrete-time 
position  profile  can  then  be  derived  for  the  target  and  fed  into  the  weapon  flyout 
model  to  compute  the  associated  measures  of  weapon  effectiveness.  The 
methodology  is  applied  to  both  land  vehicles  traveling  on  specific  roads, 
differentiated  based  on  the  road  profiles,  and  small  boats  moving  on  the  open 
sea  to  generate  various  sets  of  measures  of  weapon  effectiveness  for  different 
scenarios.  A  program  was  also  developed  to  enable  a  quick  classification  of  the 
target’s  expected  route  of  advancement  into  one  of  the  six  road  types.  This 
allows  the  user  to  perform  quick  analysis  on  the  target’s  planned  route  of  travel 
and  produce  measures  of  weapon  effectiveness,  which  can  then  be  readily 
utilized  for  swift  decision  making  on  the  weapon  options  available. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


VI 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

A.  MOTIVATION . 1 

B.  BACKGROUND . 2 

C.  THESIS  OVERVIEW . 2 

D.  LITERATURE  REVIEW . 5 

1.  Road  Profile  Derivation . 5 

2.  High-Speed  Boat  Maneuver  Estimation . 8 

3.  Map  Coordinates  Transformation . 10 

II.  ROUTE  GENERATION . 15 

A.  LAND  ROUTE  GENERATION . 15 

1.  Statistical  Road  Type  Classification . 15 

2.  Route  Synthesis  by  Road  Type . 18 

B.  SEA  ROUTE  GENERATION . 25 

III.  VELOCITY  PROFILE  COMPUTATION . 27 

A.  VELOCITY  PROFILE  FOR  LAND  VEHICLE . 27 

1 .  Radius  of  Curvature  Computation . 27 

2.  Maximum  Model  Speed  Limit  Computation . 28 

B.  VELOCITY  PROFILE  FOR  A  SEA  VESSEL . 44 

IV.  TIME  DISCRETIZATION . 49 

A.  TRAVELING  TIME  COMPUTATION . 49 

B.  CONVERSION  TO  FIXED-TIME  INTERVAL . 50 

C.  NEW  DATA  GENERATION  FOR  LAND  ROUTES . 51 

D.  NEW  DATA  GENERATION  FOR  SEA  ROUTES . 55 

V.  ROUTE  DIGITIZATION . 59 

A.  SELECTING  THE  ROUTE  OF  INTEREST . 59 

B.  IMPORTING  THE  ROUTE  OF  INTEREST . 60 

VI.  ROUTE  PROCESSING  AND  ANALYSIS . 63 

A.  CONVERSION  TO  REGULAR  INTERVALS . 63 

B.  CHANGE  OF  HEADING  COMPUTATION . 66 

VII.  CONCLUSION  AND  RECOMMENDATIONS . 73 

A.  CONCLUSION . 73 

vii 


B.  RECOMMENDATIONS 


74 


APPENDIX  A.  MATLAB  CODE  FOR  SYNTHESIZED  ROUTE 

GENERATION . 75 

APPENDIX  B.  SIMULINK  MODEL  FOR  VELOCITY  COMPUTATION . 77 

APPENDIX  C.  SIMULINK  MODEL  FOR  SEA  ROUTE  GENERATION . 79 

APPENDIX  D.  MATLAB  CODE  FOR  SYNTHESIZED  ROUTE 

PROCESSING  INTO  TSPI  FORMAT . 81 

APPENDIX  E.  MATLAB  CODE  FOR  SEA  ROUTE  PROCESSING . 87 

APPENDIX  F.  MATLAB  CODE  FOR  ROUTE  ANALYSIS . 93 

APPENDIX  G.  VBA  CODE  FOR  ROUTE  ANALYSIS . 99 

LIST  OF  REFERENCES . 107 

INITIAL  DISTRIBUTION  LIST . 109 


viii 


LIST  OF  FIGURES 


Figure  1 .  Overview  of  Weapon  Assessment  Process  Flow . 3 

Figure  2.  Centripetal  Acceleration  of  Vehicle  Traveling  along  Curved 

Road.  Adapted  from  [5] . 6 

Figure  3.  WGS-84  Ellipse . 11 

Figure  4.  ECEF  Geodetic  and  ENU  Coordinate  Frame.  Source:  [3] . 12 

Figure  5.  Statistical  Analysis  of  Coast  Ridge  Road.  Source:  [1 0] . 15 

Figure  6.  Statistical  Analysis  of  Carmel  Valley  Road.  Source:  [1 0] . 1 6 

Figure  7.  Statistical  Analysis  of  Highway  1 .  Source:  [1 0] . 1 6 

Figure  8.  Statistical  Analysis  of  Highway  1 01 .  Source:  [10] . 17 

Figure  9.  Comparison  of  Carmel  Valley  Road  and  Synthesized 

Reproduction . 19 

Figure  10.  Route  Synthesis  Path  Constraint . 20 

Figure  1 1 .  Comparison  of  Digitized  Carmel  Valley  Road  with  New 

Synthesized  Road . 21 

Figure  12.  Synthesized  Road  for  Type  1 . 22 

Figure  13.  Synthesized  Road  for  Type  2 . 22 

Figure  14.  Synthesized  Road  for  Type  3 . 23 

Figure  15.  Synthesized  Road  for  Type  4 . 23 

Figure  16.  Synthesized  Road  for  Type  5 . 24 

Figure  1 7.  Synthesized  Road  for  Type  6 . 24 

Figure  18.  Sea  Vessel  Predicted  Path  of  Travel  with  Evasive  Maneuvers . 25 

Figure  19.  Radius  of  Curvature  Computation  Using  Three  Known  Points . 27 

Figure  20.  Illustration  of  Deceleration  Constraint  in  Speed  Computation . 31 

Figure  21 .  Illustration  of  Acceleration  Constraint  in  Speed  Computation . 32 


IX 


Figure  22.  Computation  of  Speed  under  All  Three  Constraints  for  Points 

1-8 . 34 

Figure  23.  Computation  of  Speed  under  All  Three  Constraints  for  Points 

7-14 . 34 

Figure  24.  Speed  Profile  for  Road  Type  1 . 36 

Figure  25.  Acceleration  Profile  for  Road  Type  1 . 36 

Figure  26.  Speed  Profile  for  Road  Type  2 . 37 

Figure  27.  Acceleration  Profile  for  Road  Type  2 . 37 

Figure  28.  Speed  Profile  for  Road  Type  3 . 38 

Figure  29.  Acceleration  Profile  for  Road  Type  3 . 38 

Figure  30.  Speed  Profile  for  Road  Type  4 . 39 

Figure  31.  Acceleration  Profile  for  Road  Type  4 . 39 

Figure  32.  Speed  Profile  for  Road  Type  5 . 40 

Figure  33.  Acceleration  Profile  for  Road  Type  5 . 40 

Figure  34.  Speed  Profile  for  Road  Type  6 . 41 

Figure  35.  Acceleration  Profile  for  Road  Type  6 . 41 

Figure  36.  Comparison  of  Velocity  Profile . 43 

Figure  37.  Sea  Vessel  Maneuverability  Study . 45 

Figure  38.  Geometrical  Relationship  for  Sea  Vessel  Turn  Angle . 46 

Figure  39.  Sea  Vessel  Rudder  Angle  Inputs . 47 

Figure  40.  Sea  Vessel  Simulated  Course  of  T ravel . 47 

Figure  41 .  Kinematic  Model  for  Time  Duration  Computation . 49 

Figure  42.  Creation  of  New  Points  at  Discrete  Time  Intervals . 50 

Figure  43.  Type  1  Route  at  Discrete  Time  Intervals . 52 

Figure  44.  Type  2  Route  at  Discrete  Time  Intervals . 52 


x 


Figure  45.  Type  3  Route  at  Discrete  Time  Intervals . 53 

Figure  46.  Type  4  Route  at  Discrete  Time  Intervals . 53 

Figure  47.  Type  5  Route  at  Discrete  Time  Intervals . 54 

Figure  48.  Type  6  Route  at  Discrete  Time  Intervals . 54 

Figure  49.  Original  Sea  Vessel  Route . 56 

Figure  50.  Sea  Vessel  Route  at  Discrete  Time  Interval . 56 

Figure  51 .  Selection  of  Route  of  Interest . 60 

Figure  52.  Generation  of  Route  Data  in  GPX  Format . 61 

Figure  53.  Creation  of  New  Points  at  Regular  Distance  Interval . 64 

Figure  54.  Original  Route  Points  Imported . 65 

Figure  55.  New  Regularly  Spaced  Route  Points . 65 

Figure  56.  Computation  of  Pleading . 66 

Figure  57.  Computation  of  Change  in  Pleading . 67 

Figure  58.  Regularly  Spaced  Points  on  Test  Route . 67 

Figure  59.  Pleading  of  Points  on  Test  Route . 68 

Figure  60.  Change  of  Pleading  of  Points  on  Test  Route . 68 

Figure  61 .  Processed  North  Korean  Route — Type  2 . 69 

Figure  62.  Processed  Syrian  Route — Type  2 . 70 

Figure  63.  Processed  Afghanistan  Route — Type  3 . 70 

Figure  64.  Route  Analysis  in  VBA  Environment . 71 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


XII 


LIST  OF  TABLES 


Table  1 .  Statistical  Analyses  of  Actual  Roads . 17 

Table  2.  Road  Type  Classification  by  Change  in  Road  Heading . 1 8 

Table  3.  Final  Road  Type  Classification . 20 

Table  4.  Types  of  Land  Vehicles  of  Interest . 30 

Table  5.  Inputs  Characteristics  of  Sea  Vessels.  Adapted  from  [13],  [14] . 44 


xiii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


XIV 


LIST  OF  ACRONYMS  AND  ABBREVIATIONS 


CEP 

Circular  Error  Probable 

ECEF 

Earth-Centered-Earth-Fixed 

ENU 

East-North-Up 

GPS 

Global  Positioning  System 

IFS 

In-flight  Simulation 

TSPI 

Time-Space-Position  Information 

URL 

Uniform  Resource  Locator 

VBA 

Visual  Basic  for  Applications 

WGS-84 

World  Geodetic  System  -  84 

xv 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


XVI 


ACKNOWLEDGMENTS 


First,  I  wish  to  thank  my  advisors,  Professor  Morris  Driels  and  Professor 
Isaac  Kaminer,  for  their  support  and  guidance  throughout  this  project.  Their 
wealth  of  knowledge  and  sophisticated  approaches  have  been  both  very  helpful 
and  enlightening  during  the  course  of  my  thesis  research. 

In  addition,  I  would  like  to  express  my  gratitude  to  all  my  class  instructors 
and  professors  who  have  helped  me  in  one  way  or  another  during  the  course  of 
my  studies  at  NPS.  They  have  provided  me  with  an  enriching  and  stimulating 
experience  with  their  vast  knowledge  and  continual  encouragement  during  this 
long  journey. 

I  would  also  like  to  thank  my  fellow  schoolmates  at  NPS  who  showed  a 
generous  amount  of  hospitality  to  help  my  family  to  settle  seamlessly  into  a  new 
environment.  The  friendships  forged  during  my  stay  here  have  certainly  helped 
me  while  I  have  been  far  away  from  my  home  country. 

Last,  but  certainly  not  least,  I  wish  to  express  my  gratitude  toward  my 
family  for  their  utmost  support  during  the  course  of  my  study,  without  which  I 
would  not  have  been  able  to  stay  the  course. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


I.  INTRODUCTION 


A.  MOTIVATION 

In  modern  warfare,  with  the  rapid  development  of  military  technologies, 
armed  forces  around  the  world  are  increasingly  using  guided  weapons  to  protect 
their  valuable  assets  and  to  destroy  targets  of  interest  for  strategic  gains.  Such 
weapons  can  be  launched  from  various  platforms,  including  aerial,  surface,  or 
even  submerged  ones.  Even  with  the  proliferation  of  the  guided  weapons,  the 
ability  of  the  weapons  to  hit  moving  targets  remains  a  key  area  of  interest. 
According  to  [1],  the  U.S.  Navy  is  in  the  midst  of  equipping  the  submarine  and 
the  Tomahawk  missile  with  enhanced  seeker  technology  in  order  to  better  enable 
the  weapon  to  take  out  mobile  naval  threats  of  the  hostile  forces.  Meanwhile,  [2] 
reveals  that  even  the  F-35,  despite  being  the  state-of-the-art  fifth  generation 
stealth  fighter,  is  still  looking  for  options  to  enhance  its  capability  to  hit  moving 
targets. 

Depending  on  the  type  of  target  environment,  and  the  target’s 
vulnerabilities,  extensive  weapon  assessment  is  required  to  evaluate  the  likely 
weapon  performance  against  such  target,  and  the  weapon  delivery  profile  for  the 
best  possible  outcome.  To  perform  such  assessment  with  models  of  higher 
fidelity,  it  is  important  to  achieve  a  realistic  representation  of  the  target  motion 
based  on  the  target’s  type  of  environment. 

This  thesis  focuses  on  the  position  estimation  for  both  a  fast-moving  land 
vehicle  along  an  expected  route  and  a  fast-moving  boat  on  an  open  sea.  The  aim 
of  the  thesis  is  to  develop  a  methodology  to  model  an  estimator  for  the  likely  path 
of  travel  for  the  land  and  sea  targets  of  interest,  and  for  the  land  target,  to  classify 
the  target  environment  based  on  the  road  profile  and  vehicle  characteristics; 
using  this  information  enables  a  set  of  weapon  performance  estimates  to  be 
computed  to  facilitate  quick  decision  making  at  the  operational  level.  The  goal  of 
this  thesis  is  to  develop  a  realistic  model  characterizing  the  most  probable  paths 
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of  interest  for  specific  target  types  based  on  the  target  environment  to  provide 
useful  inputs  that  can  be  used  to  generate  weapon  performance  metrics  that 
enable  the  choice  of  weapon  delivery  for  best  possible  outcome. 

B.  BACKGROUND 

Weaponeering  is  defined  in  [3]  as  the  process  of  “determining  the  quantity 
of  a  specific  type  of  weapon  required  to  achieve  a  defined  level  of  target  damage, 
considering  target  vulnerability,  weapon  effects,  munitions  delivery  error,  damage 
criteria,  probability  of  kill,  weapon  reliability,  etc.”  By  determining  the  statistical 
estimate  of  the  weapon’s  effectiveness  against  specific  target,  the  user  can  then 
designate  the  appropriate  type  and  quantity  of  weapons  to  be  used  to  achieve 
the  level  of  damage  required.  For  weaponeering  to  achieve  a  close  estimate  of 
the  expected  weapon  performance,  it  is  important  for  the  parameters  assumed  to 
match  those  used  during  weapon  deployment  in  combat.  Otherwise,  an 
inaccurate  model  will  lead  to  a  poor  weaponeering  assessment,  and  result  in 
inefficient  or  ineffective  use  of  the  weapons. 

The  area  of  interest  is  in  making  realistic  weaponeering  assessment 
against  moving  targets.  Broadly,  this  can  be  classified  into  two  scenario  types, 
constrained  and  unconstrained.  The  first  scenario  type  involves  targets  moving 
along  a  defined  path,  for  example,  a  vehicle  going  down  a  particular  stretch  of 
road  or  a  train  on  a  railroad  track,  while  the  latter  scenario  involves  targets 
moving  in  an  open  environment,  for  example,  a  tank  traversing  a  desert  terrain  or 
a  speedboat  traveling  in  an  open  sea. 

C.  THESIS  OVERVIEW 

For  this  thesis,  the  focus  is  on  both  scenario  type,  particularly 
weaponeering  assessment  against  specific  type  of  land  vehicles  moving  along  a 
certain  stretch  of  road  and  small  fast-moving  boats  in  the  open  sea.  An  overview 
of  the  weapon  assessment  process  is  illustrated  in  Figure  1 . 
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Figure  1 .  Overview  of  Weapon  Assessment  Process  Flow 

Ideally,  to  perform  a  realistic  weapon  assessment,  the  ability  to  precisely 
predict  the  exact  position  of  the  target  of  specific  vehicle  type  will  enable  the 
most  accurate  simulation  of  the  actual  combat  scenario,  and  thereby  generate 
the  weapon  metrics  to  best  represent  the  likely  outcome  with  the  highest  fidelity. 
Nevertheless,  given  the  large  amount  of  time  required  for  the  weapon  In-Flight 
Simulation  (IFS)  model  to  produce  the  simulation  results  for  each  and  every  type 
of  weapon  for  a  particular  mission,  it  would  be  impractical  to  have  this  carried  out 
during  the  actual  mission  phase  when  time  is  limited. 

Therefore,  for  the  case  of  the  land  target,  it  is  necessary  to  categorize  the 
wide  distribution  of  roads  into  discrete  road  types  from  which  a  random  road, 
which  is  representative  of  each  of  the  road  types,  is  used  in  the  weapon  IFS 
model  for  the  full  weapon  assessment  to  be  carried  out.  Similarly,  for  the  case  of 
the  sea  target  in  an  unconstrained  environment,  it  is  necessary  to  predict  the 
likely  evasive  maneuvers  that  the  boat  will  adopt  in  its  approach,  and  carry  out 
the  full  weapon  assessment  in  the  weapon  IFS  model.  This  forms  the  first  part  of 
the  thesis. 
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Due  to  the  vast  number  of  roads  available,  it  is  nearly  impossible  to  carry 
out  the  full  weapon  assessment  whenever  a  new  target  emerges.  As  such,  a 
generalization  of  the  types  of  roads  available  is  required,  and  a  representative 
road  of  each  type  is  to  be  generated  and  analyzed.  In  Chapter  II  a  discussion  of 
the  generation  of  synthesized  roads  of  different  types  based  on  type  of  road 
profile  expected  for  each  type  is  offered  for  the  land  target.  The  discussion 
addresses  how  the  various  roads  are  categorized  and  how  the  representative 
road  for  each  category  is  generated.  For  the  sea  target,  the  route  that  a  boat  is 
likely  to  follow  in  its  attempt  to  evade  weapons  during  its  escape  or  attack  is 
designed  and  generated  using  an  existing  realistic  boat  simulation  model. 

In  Chapter  III,  we  examine  how  the  velocity  profile  of  the  land  vehicle 
traveling  along  each  of  the  synthesized  roads  generated  can  be  computed  based 
on  the  road  profile.  It  is  important  to  model  the  position  and  velocity  profile  of  the 
target  vehicle  as  closely  as  possible  to  its  actual  motion  to  ensure  that  the 
simulation  results  from  the  IFS  model  will  be  nearly  identical  to  the  actual 
weapon  performance  under  actual  combat  scenario.  For  the  sea  target,  the 
velocity  profile  is  generated  as  part  of  the  simulation.  Further  processing  of  the 
data  generated  is  required,  however,  to  transform  the  data  points  into  the  Time 
Space  and  Position  Information  (TSPI)  format  necessary  for  the  weapon  IFS 
model. 

Thereafter,  Chapter  IV  contains  a  study  of  how  the  position  of  both  the 
land  and  sea  vehicles,  based  on  the  velocity  profile  computed  in  Chapter  III,  can 
be  determined  and  generated  at  discrete  time  intervals,  so  that  the  TSPI  can  be 
fed  into  the  weapon  IFS  model  for  the  weapon  simulations  to  be  carried  out.  The 
IFS  model  requires  the  TSPI  inputs  for  the  target  to  be  at  one-second  intervals, 
along  with  the  velocity  profile  during  each  time-step. 

The  weapon  IFS  model  simulations  will  not  be  part  of  this  thesis  due  to  the 
classified  database  from  which  it  draws  the  weapon  information.  The 
weaponeering  community  will  make  use  of  the  TSPI  inputs  to  generate  the 
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associated  simulation  results.  Similarly,  the  simulation  results  are  classified  and 
will  not  be  disclosed  in  this  thesis. 

The  second  part  of  the  thesis  then  focuses  on  enabling  the  users  to  make 
use  of  simulation  results  generated  earlier  for  making  weaponeering  assessment 
during  the  actual  mission  planning.  For  this  part  of  the  thesis,  the  focus  is  on  the 
land  target  scenario,  due  to  the  unconstrained  nature  of  the  sea  target  scenario. 

Chapter  V  presents  the  digitization  process  of  an  actual  route  of  interest 
that  a  target  is  likely  to  travel  on.  It  discusses  how  the  map  selection  process  can 
be  carried  out,  and  how  the  route  information  can  be  imported  into  a  program 
and  processed  into  a  digital  form  that  can  be  used  for  further  target  velocity 
profile  analysis. 

Chapter  VI  then  makes  use  of  the  digitized  route  to  compute  the  likely 
velocity  profile  of  the  target,  using  such  information  to  determine  the  route  profile 
and  classify  the  route  of  interest  into  one  of  the  road  type  classification  discussed 
in  Chapter  II.  By  classifying  the  route  into  one  of  the  road  types,  the  associated 
weapon  metrics  for  the  particular  road  type  can  then  be  referred  to  for  a  close 
estimate  of  the  weapon  metrics  obtained  from  the  elaborate  weapon  simulations 
already  carried  out  for  the  specific  road  type. 

Chapter  VII  summarizes  the  conclusions  from  this  study  and  offers 
recommendations  for  further  studies. 

D.  LITERATURE  REVIEW 

In  this  section,  we  discuss  the  results  of  the  review  carried  out  on  existing 
literature  for  the  following  areas  to  aid  the  study. 

1.  Road  Profile  Derivation 

To  compute  a  realistic  velocity  estimation  for  a  vehicle  traveling  down  a 
particular  stretch  of  road,  it  is  necessary  to  relate  the  maximum  design  speed 
limit  with  the  road  curvature  along  various  stretches  of  the  road.  And  to  do  that,  it 

is  necessary  to  compute  the  curvature  at  various  points  along  the  entire  route. 
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Andersson  and  Aronsson  [4]  discuss  the  road  shape  modeling  concept  and  the 
rationale  for  the  selection  of  the  spline  interpolation  of  the  points.  The  splines 
constructed  are  approximations,  however,  and  the  splines  may  not  completely 
follow  the  exact  path  of  interest  despite  the  degree  of  shaping  done.  Hence,  an 
alternate  method  of  determining  the  road  profile  is  required. 

With  the  digitization  of  the  route,  the  route  will  be  well  defined  by  many 
points  along  the  entire  stretch  of  the  path.  Thus,  an  exact  method  of  finding  the 
curvature  at  each  and  every  point  along  the  path  can  be  used.  By  evaluating  the 
circle  that  passes  through  the  point  and  its  immediate  two  adjacent  points,  it  is 
possible  to  determine  the  curvature  of  the  road.  This  will  be  discussed  further  in 
Chapter  III. 

Once  the  curvature  of  the  path  has  been  obtained,  the  speed  limit  at  each 
point  along  the  road  can  then  be  computed.  Donnell  et  al.  [5]  discuss  the 
relationship  between  the  horizontal  road  curvature  and  the  superelevation.  By 
considering  the  centripetal  acceleration  of  the  vehicle  when  it  travels  along  a 
curved  road,  as  shown  in  Figure  2,  the  velocity  of  the  vehicle,  as  a  function  of  the 
road  curvature,  the  side  friction  factor  and  the  road  superelevation,  can  be 
determined. 


Figure  2.  Centripetal  Acceleration  of  Vehicle  Traveling  along  Curved 

Road.  Adapted  from  [5]. 
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Resolving  forces  to  both  the  horizontal  and  vertical  directions,  and  equating 
forces  in  the  horizontal  direction  only,  we  have 


Fcos0+  Nsin0  = - .  (1.1) 

R 

Similarly,  for  the  vertical  direction, 

N  cos  0  -  W  -  F  sin  0  =  0  .  (1.2) 

The  relation  between  F  and  R  is  given  by 

F  =  juN  .  (1.3) 


Using  Equation  (1.3),  and  equating  R  in  both  (1.1)  and  (1.2), 


w 


mV 2 


cos  6-  jU  sin  0  R(jU  cos  0  +  sin  0) 


Solving  for  V, 


1  -  n  tan  0 

With  Equation  (1.5),  the  following  equation  is  obtained: 


V 

15  R 


=  0.0  le  +  f. 


where 


(1.4) 


(1.5) 


(1.6) 


e  =  rate  of  superelevation,  in  percent 
f  =  side  friction  (demand  factor) 

V  =  vehicle  speed,  in  mph 
R  =  radius  of  curve,  in  feet 

With  Equation  (1 .6),  the  design  speed,  limited  by  the  road  curvature,  can 
thus  be  determined  and  incorporated  as  one  of  the  constraints  on  the  vehicle  in 
subsequent  velocity  profile  analysis. 


7 


2.  High-Speed  Boat  Maneuver  Estimation 

In  order  to  more  accurately  predict  the  evasive  maneuvers  that  the  boat 
will  adopt  in  its  approach,  it  will  be  essential  to  first  determine  how  maneuverable 
the  vessel  is.  The  maneuverability  will  ultimately  determine  the  type  of  turns  the 
vessel  is  capable  of,  as  well  as  the  time  required  and  associated  speed 
attainable  with  such  turns. 


Aslan  [6]  discusses  the  maneuvering  characteristics  of  some  of  the  fast 
craft  employed  in  different  parts  of  the  world.  Incorporating  the  correlation 
derived  by  Lewandowski  [7]  between  the  fast  craft  turning  radius  with  the  length 
LWL,  volume  V,  volumetric  Froude  number  Fv  and  rudder  angle  8,  based  on  the 
empirical  data  from  the  works  of  Denny  and  Hubble  [8],  the  following  equation 
was  developed. 


K 

LWL 


~  0.5 


1.7  +  0.0222FV 


f  \ 

2.85 

f30l 

LWL 

1 

J 

UJ 

(1.7) 


The  Froude  number,  in  turn,  is  derived  from  the  vessel  approach  speed  Ua 
and  the  volume  V  through  the  following  equation. 


Fv  = 


(1.8) 


The  normalized  turning  radius  ^  computed  and  the  approach  speed  Ua 

can  then  be  used  to  derive  the  speed  during  the  turn  Uc  through  the  following 
equation  by  Denny  and  Hubble  [8]. 


U„  =U. 


1 +  K 


/  \  2 

'  LWL  ' 


V  Rc  j 


(1.9) 


The  empirical  constant  Kc  can  be  obtained  using  the  Froude  number  Fv 
and  rudder  angle  8  using  the  following  equation. 
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(1.10) 


To  ensure  the  stability  of  the  vessel  during  a  turn,  the  Nomoto  K-T  model 
is  studied  in  [6].  To  simplify  the  model,  the  first  order  transfer  function  by  Nomoto 
was  adopted.  When  expressed  in  time  domain,  the  following  equation  is  derived 
and  used  in  the  model. 


where 


Tr  +  r  =  KS . 


r  =  lf/ . 


(1.11) 

(1.12) 


Assuming  the  yaw  inertia  coefficient  T  to  be  zero,  Aslan  [6]  then  derives 
the  following  relation  between  the  turning  ability  K  and  the  vessel  approach 
speed  Ua,  the  turning  radius  Rc  and  the  max  rudder  angle  Smax. 


K  = 


U„ 


RxS,„ 


(1.13) 


Using  the  computed  value  of  K,  substituting  Equation  (1.13)  into  Equation 
(1.1 1)  and  introducing  a  correction  factor  ij  of  value  0.909  to  ensure  consistency 
with  the  turning  radius  Rc  in  the  Lewandowski’s  equation  [7],  the  following 
equation  is  obtained. 

\j/  =  Kdri.  (1-14) 


Equation  (1.14)  can  then  be  integrated  with  respect  to  time  for  the 
computation  of  the  velocity  profile  using  the  following  equations. 


x  =  Uc  sin  Y 
y  =  Uc  cos  y 


(1.15) 


Further  integration  of  Equation  (1.15)  can  generate  the  position  profile  of 
the  vessel  with  respect  to  time. 
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3.  Map  Coordinates  Transformation 

To  digitize  a  road,  it  is  necessary  to  convert  map  locations  in  the  Geodetic 
Coordinates  System  to  that  of  a  local  reference  frame.  Both  Farrell  [9]  and 
Driels  [3]  discuss  the  conversion  of  the  Global  Positioning  System  (GPS) 
readings  from  the  WGS-84  geodetic  system  to  a  local  East-North-Up  (ENU) 
Cartesian  axis. 

The  basic  reference  frame  for  terrestrial  navigation  is  the  Earth-Centered- 
Earth-Fixed  (ECEF)  geodetic  coordinate  system,  where  the  coordinates  are  in 
latitude  0,  longitude  A,  and  altitude  h.  This  is  also  referred  to  as  the  global 
reference  frame,  which  is  most  commonly  associated  with  the  GPS. 

An  alternate  reference  frame  is  known  as  the  ECEF  rectangular  system, 
where  the  x-axis  passes  from  the  origin  at  the  Earth’s  center  to  the  prime 
meridian,  the  z-axis  passes  through  the  North  Pole,  and  the  y-axis  orthogonal  to 
these  two  axes.  This  frame  is  convenient,  as  it  allows  for  easy  relation  to  the 
ENU  frame  that  is  used  subsequently. 

To  use  the  ECEF  geodetic  system,  the  shape  of  the  Earth’s  ellipsoid 
needs  to  be  defined,  and  the  parameters  for  the  assumed  shape  used  in  GPS, 
known  as  the  World  Geodetic  System  -  1984  (WGS-84),  are  defined  by  the 
equatorial  and  polar  radii,  a  and  b,  of  6,3778,137.0  m  and  6,356752.3142  m, 
respectively.  The  ellipse  is  presented  in  Figure  3. 
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Figure  3.  WGS-84  Ellipse 


The  ellipsoid  flatness  is  defined  by 

/=— ■  (1.16) 

a 

The  eccentricity  of  the  ellipsoid  is  defined  by 

*  =  V/(2-/)-  (1-17) 

The  distance  from  the  surface  of  the  ellipsoid  to  the  z-axis  along  the 
ellipsoid  normal,  as  a  function  of  0,  is  defined  by 

W)=  1  =.  2  ■  0-18) 

y  1-e  sin  (j) 

The  ECEF  rectangular  coordinates  can  then  be  calculated  from  the 
geodetic  system  using  the  following: 


x  =  (h  +  N)cost/)cosA, 

(1.19) 

y  =  (h  +  N)cos0sinA  , 

(1.20) 
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z  =  (/z  +  (l-e2)iVjsin0. 


(1.21) 


Thereafter,  in  order  to  form  a  right-handed  coordinate  system  that  is 
analogous  to  the  Cartesian  coordinate  system,  it  is  essential  to  further  transform 
from  the  ECEF  rectangular  system  to  that  of  the  ENU  system.  This  rectangular 
system  is  defined  with  the  origin  at  a  point  of  interest  on  the  Earth’s  surface,  with 
the  vertical  axis  along  the  ellipsoid  normal,  and  the  other  two  axes  orthogonal 
and  aligned  to  the  local  geographic  East  and  North  directions.  The  ENU  is 
presented  in  Figure  4. 


Tecef 


Figure  4.  ECEF  Geodetic  and  ENU  Coordinate  Frame.  Source:  [3]. 

With  the  origin  of  the  ENU  defined  as  (x0,y0,z0)  in  the  ECEF  rectangular 
frame,  to  transform  the  ECEF  rectangular  frame  to  the  ENU  frame  would  require 
a  translation  of  the  origin  to  (x0,y0,z0),  followed  by  a  rotation  of  the  axes  to  align 
with  those  of  the  ENU.  The  complete  transformation  is  given  by 
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(1.22) 


-sin  A 

cos  A 

0 

[x-Xo 

n 

- 

-  sin  (f>  cos  A 

-  sin  (f>  sin  A 

COS0 

V! 

1 

V! 

o 

Ku) 

cos  (f)  cos  A 

cos  (f)  sin  A 

sin^ 

l - 

1 

o 

With  the  transformation  in  Equation  (1.22)  into  the  ENU  system,  which  is 
analogous  to  the  familiar  Cartesian  system,  the  digitized  route  information  can 
then  be  easily  further  processed  in  subsequent  analysis. 
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II.  ROUTE  GENERATION 


A.  LAND  ROUTE  GENERATION 

In  this  section,  the  aim  is  to  generate  routes  that  are  representative  of  the 
different  road  types.  Before  we  begin  to  generate  synthesized  road  for  each  road 
type,  it  is  essential  to  first  define  a  quantitative  approach  to  classify  the  road 
types.  Thereafter,  we  will  generate  the  roads  based  on  the  characteristic 
parameter  associated  with  each  road  type. 

1.  Statistical  Road  Type  Classification 

Driels  [10]  introduced  the  concept  of  road  type  characterization  through 
the  analysis  of  actual  roads.  A  statistical  analysis  is  carried  out  to  study  both  the 
heading  and  change  in  heading  for  four  distinct  types  of  real  roads  to  determine 
their  respective  means  and  standard  deviation.  The  statistical  results,  in 
histogram  form,  are  presented  in  Figures  5  to  8. 


350  -| 


Coast  Ridge  Road 


120.00% 


Heading  Mean  (deg)  121.036 
Heading  StDev  (deg)  54.639 


300  -  A  Heading  Mean  (deg)  -0.020 
A  Heading  StDev  (deg)  13.419 
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Figure  5.  Statistical  Analysis  of  Coast  Ridge  Road.  Source:  [10]. 
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Figure  6.  Statistical  Analysis  of  Carmel  Valley  Road.  Source:  [10]. 
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Figure  7.  Statistical  Analysis  of  Highway  1.  Source:  [10]. 
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Hwy  101 


Heading  Mean  (deg)  138.583 
Heading  StDev  (deg)  22.086 
A  Heading  Mean  (deg)  -0.009 
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Figure  8.  Statistical  Analysis  of  Highway  101.  Source:  [10]. 

The  results  of  the  preceding  statistical  analyses  are  consolidated  in  Table 
1  for  ease  of  comparison. 


Table  1 .  Statistical  Analyses  of  Actual  Roads 


Road  type 

Mean 

heading 

Sigma 

heading 

Mean 
change  in 
heading 

Sigma 
change  in 
heading 

Coast  Ridge 
Road 

121 

54.6 

-0.02 

13.42 

Carmel  Valley 
Road 

120 

43.3 

-0.01 

4.93 

Pacific  Coast 
Hwy  (Hwy  1 ) 

138 

31.0 

-0.10 

3.63 

Highway  101 

139 

21.1 

-0.01 

0.49 
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From  the  results  obtained,  it  can  be  easily  observed  that  the  change  in 
road  heading  generally  follows  a  normal  distribution,  and  the  standard  deviation 
in  the  change  of  heading  increases  with  the  number  of  curves  and  degree  of 
curvatures  expected  on  the  roads.  From  this,  for  the  synthesized  road 
generation,  the  synthesized  roads  are  assumed  to  follow  a  normal  distribution 
and  parameterized  by  the  standard  deviation  of  the  road  change  in  heading.  A 
total  of  six  different  road  types  are  identified,  and  their  associated  standard 
deviations  in  change  of  heading  are  summarized  in  Table  2. 


Table  2.  Road  Type  Classification  by  Change  in  Road  Fleading 


Road 

type 

Sigma  change  in 
heading 

Mean  sigma  change  in 
heading 

1 

0<  to  <1 

0.5 

2 

1<  to  <4.5 

2.5 

3 

4.5<  to  <9.5 

7 

4 

9.5<  to  <14.5 

12 

5 

14. 5<  to  <19.5 

17 

6 

19.5<  to  <30.5 

25 

2.  Route  Synthesis  by  Road  Type 

Using  the  parameters  defined  in  Table  2  alone  is  not  enough  to  fully 
characterize  and  produce  synthesized  roads  that  are  representative  of  the 
original  road.  A  comparison  of  the  original  Carmel  Valley  and  the  synthesized 
road  reproduced  using  the  same  road  parameters  is  presented  in  Figure  9. 
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Figure  9.  Comparison  of  Carmel  Valley  Road  and 
Synthesized  Reproduction 

It  can  be  clearly  seen  that  the  synthesized  reproduction  of  the  Carmel 
Valley  road  is  totally  random  and  lacks  a  general  direction  of  travel.  Therefore, 
synthesizing  roads  in  such  an  unconstrained  manner  is  not  suitable,  and 
additional  constraints  are  required  to  produce  synthesized  roads  that  can  be 
more  representative  of  the  original.  To  prevent  excessive  randomness  in  the  path 
generation  and  to  impose  a  general  road  heading,  the  mean  heading  and  the 
standard  deviation  in  the  heading  are  included  for  the  road  synthesis. 

The  road  synthesis  is  carried  out  by  iteratively  generating  small  segments 
of  path  lengths  with  change  of  heading  randomly  drawn  from  a  normal 
distribution  with  the  mean  and  parameters  as  defined  in  Table  3.  These  road 
segments  are  then  incrementally  added  on  to  an  arbitrarily  chosen  starting  point 
until  a  finite  route  distance  is  achieved.  In  addition  to  this,  a  further  constraint  is 
also  imposed  on  the  road  heading  by  restricting  the  maximum  change  in  heading 
to  no  more  than  two  standard  deviations  of  the  heading.  This  is  illustrated  in 
Figure  10. 
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Table  3.  Final  Road  Type  Classification 


Road 

type 

Mean  heading 

Sigma 

heading 

Sigma  change  in 
heading 

Mean  sigma  change 
in  heading 

1 

135 

25 

0<  to  <1 

0.5 

2 

135 

29 

1<  to  <4.5 

2.5 

3 

135 

41 

4.5<  to  <9.5 

7 

4 

135 

53 

9.5<  to  <14.5 

12 

5 

135 

66 

14. 5<  to  <19.5 

17 

6 

135 

73 

19.5<  to  <30.5 

25 

With  these  additional  constraints  imbued  into  the  program,  a  new 
synthesized  road  is  generated  for  verification  of  the  updated  program.  The  new 
synthesized  road,  while  not  exactly  the  same  as  the  actual  Carmel  Valley  road, 
which  is  expected  due  to  the  randomness  in  the  generation  of  paths  during  the 
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road  synthesis,  is  much  closer  to  the  actual  road  than  before,  and  more 
representative  of  a  road  in  the  specific  road  type.  The  roads  are  provided  in 
Figure  11  for  comparison. 


Figure  1 1 .  Comparison  of  Digitized  Carmel  Valley  Road  with  New 

Synthesized  Road 


With  the  new  changes  incorporated,  and  having  verified  the  ability  of  the 
program  to  produce  a  synthesized  road  that  is  reasonably  representative  of 
actual  digitized  road,  the  representative  synthesized  roads  for  each  of  the  six 
distinct  road  types  are  generated  and  presented  in  Figures  12  to  17.  The 
Cartesian  coordinates  of  all  the  points  making  up  the  newly  synthesized  roads 
are  saved  for  subsequent  work  in  Chapter  III  for  conversion  into  the  TSPI  input 
format  required  for  weapon  performance  analysis  by  the  IFS  model. 

As  intended,  the  roads  generated  all  have  a  common  heading  towards  the 
southwestern  direction,  and  the  road  profile  gets  increasingly  winding  as  the 
standard  deviation  in  the  change  of  heading  increases.  From  the  plots  of  the 
routes  generated,  it  can  be  easily  observed  that  the  routes  are  as  per  design, 
and  they  are  ready  to  be  worked  on  further  in  subsequent  chapters.  A  copy  of  the 
MATLAB  codes  for  the  generation  of  routes  for  the  various  road  types  can  be 
found  in  Appendix  A. 
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Type  1 


Easting  (km) 

Figure  12.  Synthesized  Road  for  Type  1 


Figure  13.  Synthesized  Road  for  Type  2 
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Northing  (km)  Northing  (km) 


Type  3 


Figure  14.  Synthesized  Road  for  Type  3 


Type  4 


Figure  15.  Synthesized  Road  for  Type  4 
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Northing  (km)  Northing  (km) 


Type  5 


Figure  1 6.  Synthesized  Road  for  Type  5 


Type  6 


Figure  1 7.  Synthesized  Road  for  Type  6 
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B.  SEA  ROUTE  GENERATION 


In  earlier  discussion,  the  target’s  route  of  travel  is  constrained  by  a  specific 
path,  and  these  paths  can  in  turn  be  characterized  by  their  profile  into  different 
road  types,  from  which  a  representative  path  can  be  generated  for  further 
analysis.  Yet,  for  a  vessel  moving  in  the  open  sea,  there  is  no  such  constraint  to 
restrict  it  within  a  fixed  path.  Instead,  it  is  free  to  move  wherever  it  wants, 
whenever  it  wants.  In  order  to  generate  a  realistic  prediction  of  the  target’s  path 
of  travel  on  the  open  sea,  certain  assumptions  about  its  travel  are  required. 

For  this  study,  it  is  assumed  that  the  general  heading  for  the  target,  when 
viewed  on  an  X-Y  Cartesian  plane,  is  towards  the  right,  along  the  x-axis.  To 
reach  its  target  in  the  shortest  amount  of  time,  it  always  maintains  its  approach 
speed  at  the  maximum  possible.  To  try  to  evade  a  weapon  directed  at  it,  though, 
the  vessel  must  execute  turns  at  the  highest  turning  speed  possible,  while 
maintaining  its  general  motion  towards  the  desired  heading.  This  is  termed  as  the 
desired  travel  path  for  the  sea  vessel,  as  illustrated  in  Figure  18. 


Figure  18.  Sea  Vessel  Predicted  Path  of  Travel  with  Evasive  Maneuvers 
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Of  course,  a  vessel  cannot  perform  sharp  turns,  as  depicted  by  the 
desired  travel  path  in  Figure  18,  due  to  the  turning  radius  required  when  making 
turns  at  high  speeds.  As  such,  a  more  realistic  path  for  a  sea  vessel  executing 
evasive  maneuvers  at  high  speeds  is  more  accurately  represented  by  the  blue 
path,  as  shown  in  Figure  18. 

Furthermore,  for  the  turning  angle  0,  the  sharper  the  angle,  the  higher  is 
the  speed  of  the  sea  vessel  along  the  x-direction.  Thus,  a  large  angle  will  enable 
the  target  to  move  down  the  x-axis  in  the  shortest  possible  time  while  executing 
the  evasive  maneuvers.  A  large  angle  turn,  however,  will  not  be  aggressive 
enough  for  an  evasive  maneuver.  Hence  a  sharp  turn  angle  of  60  degrees  is 
used  in  the  predicted  route  of  the  sea  vessel  to  simulate  more  aggressive 
evasive  maneuver  while  maintaining  a  substantial  speed  along  the  x-axis  for  the 
purpose  of  this  study. 
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III.  VELOCITY  PROFILE  COMPUTATION 


A.  VELOCITY  PROFILE  FOR  LAND  VEHICLE 

In  this  section,  the  aim  is  to  determine  the  velocity  profile  of  the  land 
vehicle  traveling  along  a  constrained  path.  The  model  maximum  speed  of  the 
vehicle  is  constrained  by  both  the  vehicle  maximum  speed  and  the  road  design 
speed.  As  mentioned  in  Chapter  I,  in  order  to  make  use  of  Equation  (1.6)  to 
compute  the  maximum  road  design  speed  limit  along  a  curve,  it  is  necessary  to 
determine  the  radius  of  curvature. 

1 .  Radius  of  Curvature  Computation 

Six  distinct  routes,  made  up  of  thousands  of  intermediate  points  with  their 
associated  Cartesian  coordinates,  have  been  generated  earlier.  Due  to  the  large 
number  of  points  in  each  route,  it  is  assessed  that  it  will  be  computationally  less 
intensive  to  solve  for  the  radius  of  curvature  at  a  point  using  a  circle  that  passes 
through  the  point  and  the  two  points  immediately  adjacent  to  it,  as  shown  in 
Figure  19.  Furthermore,  this  method  produces  an  exact  solution  for  the  radius  of 
curvature,  while  there  may  be  certain  loss  of  accuracy,  even  with  good  shaping, 
with  the  method  of  spline  interpolation,  as  suggested  by  Andersson  [4]. 


(x2,y2) 


Figure  19.  Radius  of  Curvature  Computation  Using  Three  Known  Points 
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The  gradients  of  the  lines  connecting  the  first  and  second  pair  of  points 
are  given  by 


X2  -X, 


x3  -x2 


(3.1) 


The  equation  of  the  lines  connecting  the  midpoint  of  line  of  gradient  ma 
and  that  of  gradient  mb  to  the  center  of  the  circle  of  curvature  are  given  by 
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Equating  the  value  of  yc  using  the  equation  of  lines  in  (3.2),  we  get 
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Using  Equations  (3.2)  and  (3.3),  the  center  of  the  circle  of  curvature,  and 
hence  the  radius  of  curvature,  can  be  computed  for  every  point  on  the 
synthesized  route,  less  the  start  and  end  points.  Having  determined  the  radius  of 
curvature,  we  can  then  commence  on  establishing  the  velocity  profile  of  a  target 
traveling  along  the  designated  route. 


2.  Maximum  Model  Speed  Limit  Computation 

Having  computed  the  radius  of  curvature  in  the  previous  section,  we  can 
then  proceed  to  compute  the  maximum  road  design  speed  at  each  point  of  the 
route  based  on  the  road  profile  using  the  Equation  (1 .6)  discussed  by  Donnell  [5], 
as  mentioned  in  Chapter  II.  Ideally,  the  vehicle  should  travel  at  the  road  design 
speed  limit  in  order  to  minimize  the  time  taken  to  reach  the  destination.  However, 
this  is  only  possible  if  either  the  route  that  the  vehicle  is  traveling  on  is  fairly 
straight  with  little  change  on  the  maximum  speed  limit,  or  the  vehicle  has 
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extraordinary  acceleration  and  braking  capabilities  to  enable  it  to  change  to  the 
desired  speed  at  will.  Even  a  Formula  1  race  driver  racing  down  a  straight  stretch 
of  road  must  brake  hard  and  early  to  successfully  negotiate  a  sharp  bend  ahead 
without  skidding  off  the  road,  and  thereafter  accelerate  hard  down  the  straight 
road  to  reach  maximum  speed.  Therefore,  to  calculate  the  realistic  maximum 
model  speed  that  the  vehicle  is  able  to  travel  along  any  route,  not  only  is  it 
important  to  determine  the  road  design  speed  limit  to  prevent  the  vehicle  from 
skidding  out  during  a  turn,  it  is  also  equally  important  to  take  into  consideration 
the  vehicle’s  performance  capabilities,  i.e.,  the  vehicle’s  maximum  speed,  as  well 
as  its  acceleration  and  deceleration  limits. 

To  determine  the  maximum  model  speed  for  the  vehicle  at  any  point  along 
the  route  of  travel,  it  is  necessary  to  compare  the  design  road  speed  limit  at  each 
point  with  the  maximum  vehicle  speed  limit,  and  choose  the  lower  of  the  two. 
This  will  set  the  maximum  model  speed  that  the  vehicle  cannot  exceed  at  every 
point  along  the  road. 

For  the  purpose  of  weaponeering  assessment,  there  are  generally  three 
classes  of  vehicles  that  are  of  interest.  For  this  thesis,  the  focus  is  on  the  first 
class  of  vehicle  where  the  maximum  vehicle  speed  is  expected  to  be  around  60 
mph.  For  the  vehicle  acceleration  limit,  the  time  taken  for  an  average  vehicle  to 
reach  100  km/h  is  around  ten  seconds,  and  assuming  a  linear  acceleration 
profile,  the  acceleration  limit  is  assumed  as  0.3g  for  our  study.  Assuming  that  the 
braking  is  twice  as  effective,  the  deceleration  is  assumed  as  0.6g.  The 
breakdown  of  the  expected  vehicle  performance  for  all  three  classes  is 
summarized  in  Table  4.  With  the  speed,  acceleration,  and  deceleration  limits  all 
established,  it  is  now  possible  to  compute  the  actual  speed  of  the  vehicle  at 
every  point  along  the  route  of  interest. 
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Table  4.  Types  of  Land  Vehicles  of  Interest 


Vehicle  type 

Max  speed 
ms'1  (mph) 

Max 

acceleration  (g) 

Max  braking  (g) 

Fast 

26.82  (60) 

0.3 

0.6 

Slow 

13.41  (30) 

0.15 

0.3 

X-slow 

6.71  (15) 

0.08 

0.16 

To  determine  the  actual  speed  of  the  vehicle  at  each  point  along  the  route, 
it  is  important  to  first  identify  the  points  near  the  current  point  where  the  model 
maximum  speed  is  lower  than  that  for  the  current  point.  Imagining  ourselves  as 
the  driver  of  the  vehicle  and  assuming  we  reach  the  start  of  the  route  of  interest 
with  a  certain  reasonable  starting  speed,  at  the  current  instance,  we  are  only 
concerned  with  any  point  lying  ahead  that  requires  us  to  take  on  a  lower  speed 
than  our  current  speed,  as  these  points  will  determine  the  time  and  amount  of 
braking  needed  to  reach  the  lower  speed  required  in  time. 

An  illustration  of  the  deceleration  constraint  is  provided  in  Figure  20.  The 
vehicle  is  traveling  at  27  m/s  currently,  and  a  point  60  meters  ahead  has  a  model 
maximum  speed  of  12.4  m/s,  and  the  vehicle  is  capable  of  decelerating  at  6  m/s2. 
Assuming  a  linear  deceleration  profile,  the  relationship  between  the  distance 
traveled  s,  the  initial  velocity  u,  and  the  final  velocity  v  is  given  by 

v2=u2  +  2as.  (3.4) 

Using  Equation  (3.4),  the  vehicle  is  capable  of  reaching  12.4  m/s  from  27  m/s  in 
a  distance  of  48  meters.  Thus,  there  is  no  need  for  the  vehicle  to  reduce  speed  at 
the  current  instance,  and  it  can  continue  traveling  at  the  model  maximum  speed 
until  point  2.  If  the  model  maximum  speed  at  point  2  is  higher,  it  can  even 
accelerate  to  27.4  m/s  before  decelerating  to  12.4  m/s  at  point  7.  Once  it  reaches 
point  2,  the  vehicle  can  maintain  the  model  maximum  speed,  and  it  will  need  to 
begin  decelerating  in  order  to  meet  the  model  maximum  speed  at  point  7 
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subsequently.  For  point  3,  due  to  the  lower  model  maximum  speed  at  point  7,  the 
maximum  speed  possible  is  only  25.2  m/s,  and  the  vehicle  will  need  to 
decelerate  to  this  speed  when  it  reaches  point  3. 


Figure  20.  Illustration  of  Deceleration  Constraint  in  Speed  Computation 

Thus  far,  the  focus  of  the  example  using  Figure  20  has  been  on  the  effect 
of  the  lowest  model  maximum  speed  of  12.4  m/s  at  point  7  on  the  maximum 
possible  speed  at  all  the  points  before  point  7.  Nevertheless,  it  is  important  to 
note  that  the  effects  of  this  for  all  the  points  with  lower  model  maximum  speeds 
will  need  to  be  considered  instead.  A  model  maximum  speed  of  12.5  m/s  at  point 
6,  though  still  higher  than  that  at  point  7,  would  have  a  greater  impact  on  the 
traveling  speed  at  the  earlier  points  due  to  the  shorter  stopping  distance 
available.  Thus,  the  model,  as  it  steps  through  each  point  on  the  path,  will  need 
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to  consider  all  the  points  with  lower  model  maximum  speeds  within  a  specific 
lookahead  distance,  and  compute  the  maximum  possible  speed  for  the  vehicle  at 
the  current  point  based  on  the  deceleration  limit  and  the  stopping  distance 
available  to  each  of  these  possible  deceleration  “limiting  factors.” 

Along  the  same  line  of  reasoning,  the  application  for  the  acceleration  limit 
constraint  must  also  be  considered  as  the  vehicle  travels  from  point  to  point 
sequentially  along  the  route  of  interest.  For  the  acceleration  constraint,  the  points 
of  interest  are  those  prior  to  the  current  point  with  lower  model  maximum  speeds. 
An  illustration  of  the  acceleration  constraint  is  provided  in  Figure  21.  which  is  a 
continuation  of  the  previous  example  in  Figure  20. 


Figure  21 .  Illustration  of  Acceleration  Constraint  in  Speed  Computation 
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For  point  8,  the  maximum  speed  possible  due  to  the  acceleration  limit  is 
14.6  m/s;  thus,  the  constraint  in  this  case  is  due  to  the  model  maximum  speed. 
For  point  9,  the  maximum  speed  of  the  vehicle  at  point  7  is  only  12.4  m/s,  and 
the  vehicle  is  capable  of  accelerating  at  3  m/s2.  Assuming  a  linear  acceleration 
profile,  using  Equation  (3.4),  the  maximum  possible  speed  due  to  the 
acceleration  limit  from  point  7  over  the  20m  distance  is  16.5  m/s,  which  is  lower 
than  the  model  maximum  speed  of  17.6  m/s.  Yet,  the  realizable  speed  of  the 
vehicle  is  lower  at  16.0  m/s  instead  due  to  the  effect  of  the  model  maximum 
speed  at  point  8,  with  a  pick-up  distance  of  10  m  only. 

Thus,  the  model,  as  it  steps  through  each  point  on  the  path,  will  need  to 
consider  all  the  points  with  lower  model  maximum  speeds  within  a  specific 
lookback  distance  and  compute  the  maximum  possible  speed  for  the  vehicle  at 
the  current  point  based  on  the  acceleration  limit  and  the  pick-up  distance 
available  to  each  of  these  possible  acceleration  “limiting  factors.” 

Flaving  introduced  both  the  deceleration  and  acceleration  constraints  for 
the  model,  we  then  need  to  combine  both  these  constraints  for  the  computation 
of  the  actual  speed  of  the  vehicle.  This  is  easily  achieved  by  comparing  the 
maximum  speeds  possible  due  to  both  the  acceleration  and  deceleration 
constraints  with  the  model  maximum  speed  at  every  point  along  the  path,  and 
choosing  the  lowest  of  the  three  to  satisfy  all  the  constraints.  All  three  constraints 
are  presented  for  the  two  earlier  examples  in  Figures  22  and  23. 

From  Figure  22,  we  can  see  that  the  actual  speed  of  the  vehicle  at  points 
1-8  is  primarily  constrained  by  the  deceleration  limit  and  the  model  maximum 
speed.  This  is  expected  as  the  model  speed  limit  along  this  stretch  is  decreasing 
in  general. 
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Figure  22.  Computation  of  Speed  under  All  Three  Constraints 

for  Points  1-8 


Figure  23.  Computation  of  Speed  under  All  Three  Constraints 

for  Points  7-14 
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From  Figure  23,  we  can  clearly  see  that  the  speed  for  the  vehicle  at  points 
9  to  1 1  is  determined  by  the  acceleration  limits  of  the  vehicle.  The  computation  of 
the  speeds  at  points  10  and  12  is  of  interest,  as  the  values  of  the  speed  due  to 
the  acceleration,  deceleration,  and  model  maximum  speeds  constraints  are  all 
distinctly  different.  To  meet  the  constraints  of  all  three  limits,  the  speed  of  the 
vehicle  is  chosen  to  be  the  lowest  of  the  three  limits.  For  point  10,  the  deciding 
“limiting  factor”  is  that  of  the  acceleration  limit,  while  that  for  point  12  is  the 
deceleration  limit  determined  by  the  low  model  speed  limit  at  point  13. 

The  approach  that  we  have  discussed  thus  far  is  then  applied  to  all  the 
points  along  the  path  of  interest  to  generate  the  highest  possible  speed  of  the 
vehicle  at  each  of  the  points,  thereby  generating  the  velocity  profile  of  the  vehicle 
that  satisfies  all  three  physical  constraints  discussed  earlier. 

The  velocity  and  acceleration  profile  for  each  of  the  six  road  types  are 
plotted  along  the  entire  path  to  verify  that  the  speed  and  acceleration  constraints 
are  fully  satisfied.  Due  to  the  large  number  of  points  involved  for  each  plot,  a  plot 
of  100  sampling  points  is  provided  for  each  of  the  road  types  for  better  clarity. 
The  respective  plots  are  provided  in  Figures  24  to  35. 
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Figure  24.  Speed  Profile  for  Road  Type  1 


Figure  25. 


Acceleration  Profile  for  Road  Type  1 
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Figure  26.  Speed  Profile  for  Road  Type  2 


Figure  27.  Acceleration  Profile  for  Road  Type  2 
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Figure  28.  Speed  Profile  for  Road  Type  3 


Figure  29.  Acceleration  Profile  for  Road  Type  3 
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Figure  30.  Speed  Profile  for  Road  Type  4 


Figure  31 .  Acceleration  Profile  for  Road  Type  4 
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Figure  32.  Speed  Profile  for  Road  Type  5 


Figure  33.  Acceleration  Profile  for  Road  Type  5 
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Figure  34.  Speed  Profile  for  Road  Type  6 


Figure  35.  Acceleration  Profile  for  Road  Type  6 
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From  all  the  velocity  and  acceleration  plots  presented,  we  can  see  that  the 
speeds  of  the  vehicle  are  kept  at  less  than  or  equal  to  the  model  maximum 
speeds,  while  the  acceleration  of  the  vehicle  does  not  exceed  the  limits  imposed. 
This  implies  that  the  solution  meets  all  the  kinematic  constraints  satisfactorily;  as 
such,  the  results  obtained  from  the  approach  discussed  earlier  are  valid. 

During  the  course  of  the  thesis,  an  optimization  approach  is  explored  to 
determine  the  solution  for  a  time-minimization  problem.  Gill  [11]  discusses  the 
use  of  the  SNOPT  software  for  the  constrained  optimization  problems,  which  is 
applicable  to  this  thesis.  Numerous  optimization  runs  are  attempted  with  no 
meaningful  results  obtained  until  the  end  of  the  thesis  work.  As  further  tuning  of 
the  optimization  program  is  required  to  generate  the  results  that  meet  all  the 
constraints,  the  optimization  approach  is  brought  to  a  halt,  and  further  exploration 
of  this  area  is  recommended  for  future  work  to  obtain  the  optimal  solution  to  the 
time-minimization  problem. 

An  alternate  approach  for  the  velocity  computation  is  also  explored  toward 
the  end  of  the  thesis  to  include  the  effects  of  the  vehicle  dynamics  in  the 
computation  of  the  vehicle  velocity  profile.  Kaminer  [12]  developed  a  model  in  the 
Simulink  environment  to  include  the  effects  of  vehicle  dynamics  on  the  vehicle 
trajectory  and  its  velocity  profile  as  it  travels  down  a  particular  route. 

The  velocity  profile  computation  is  carried  out  using  the  model  developed 
by  Kaminer  [12],  and  the  results  for  the  initial  approach  and  that  obtained  using 
the  Simulink  model  are  provided  in  Figure  36.  Due  to  the  large  number  of  data 
points  involved,  a  small  sample  of  the  data  points  is  illustrated  in  the  plots  for 
ease  of  comparison. 
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Figure  36.  Comparison  of  Velocity  Profile 


From  the  obtained,  it  can  be  seen  that  the  velocity  profile  generated  by 
Simulink  does  not  produce  the  same  sharp  changes  in  velocity  that  are  observed 
in  the  velocity  profile  obtained  initially;  as  such,  the  Simulink  model  is  able  to 
better  capture  the  effects  of  vehicle  dynamics  and  produce  a  more  realistic 
representation  of  the  actual  vehicle  velocity  profile.  However,  it  can  also  be 
observed  that  the  velocity  profiles  generated  in  both  approaches  are  largely 
similar.  Therefore,  the  velocity  profile  generated  using  the  approach  discussed 
earlier  provides  a  very  good  approximation  of  the  actual  vehicle  performance 
under  the  given  constraints. 

In  view  of  the  limited  time  remaining  for  the  completion  of  the  thesis  work, 
the  Simulink  model  approach  is  not  studied  further.  As  the  results  from  the  initial 
approach  are  very  close  to  those  obtained  from  the  Simulink  model,  the  original 
set  of  velocity  profiles  computed  is  retained.  For  the  rest  of  the  thesis,  the  focus 
is  on  using  the  results  from  the  initial  approach  used  for  the  computation  of  the 
vehicle  velocity  profile.  For  future  work  in  similar  areas,  the  Simulink  model 
approach  is  strongly  recommended.  A  copy  of  the  Simulink  model  for  velocity 
computation  is  provided  in  Appendix  B. 
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B.  VELOCITY  PROFILE  FOR  A  SEA  VESSEL 


For  the  case  of  the  sea  vessel,  the  approach  is  different,  as  the  simulation 
model  by  Aslan  [6]  is  able  to  generate  the  velocity  profile  of  the  sea  vessel 
throughout  the  course  of  its  motion.  The  course  of  the  sea  vessel  is  determined 
by  the  rudder  inputs,  and  its  maneuverability  dependent  on  the  sea  vessel 
characteristics.  The  sea  vessel  characteristics  that  are  required  as  inputs  to  the 
simulation  model  for  two  examples  of  sea  vessels  of  interest  are  summarized  in 
Table  5.  Therefore,  the  aim  is  to  generate  a  suitable  set  of  rudder  input  to 
produce  the  desired  path  of  travel,  shown  in  Figure  18,  for  the  sea  vessel  and 
generate  the  associated  velocity  profile. 


Table  5.  Inputs  Characteristics  of  Sea  Vessels.  Adapted  from  [13],  [14]. 


Sea  Vessel 

Length  on 
Waterline 
(m) 

Displacement 

(ton) 

Maximum 

Rudder 

Angle 

(deg) 

Approach 

Speed 

(knots) 

Bladerunner  51 

15 

16 

30 

65 

x  * 

Boston  Whaler 


To  generate  the  correct  rudder  inputs  to  produce  the  desired  path  of  travel 
for  the  sea  vessel,  it  is  necessary  to  first  understand  the  maneuverability  of  the 
sea  vessel,  which  is  vessel  specific  and  dependent  on  the  vessel  length,  weight, 
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and  maximum  speed.  To  do  that,  we  first  set  the  rudder  angle  input  at  the 
maximum  rudder  angle  throughout  the  simulation  run.  This  will  simulate  the  sea 
vessel  making  the  tightest  turn  when  going  at  its  maximum  speed,  and  from  this, 
we  can  determine  the  time  taken  for  the  boat  to  make  a  complete  turn.  The 
results  of  the  maneuverability  study  are  shown  in  Figure  37. 


Sea  Vessel  Maneuverability  Study 


Figure  37.  Sea  Vessel  Maneuverability  Study 

Through  the  maneuverability  study,  we  determine  that  the  Bladerunner  51 
takes  around  41.05  seconds  to  turn  through  a  complete  round  of  tactical 
diameter  389  meters  at  its  maximum  turning  speed  of  29.8  m/s.  By  contrast,  the 
Boston  Whaler  takes  around  13.04  seconds  to  turn  through  a  complete  round  of 
tactical  diameter  68  meters  at  its  maximum  turning  speed  of  16.3  m/s.  With  this 
information,  we  can  then  proceed  to  work  on  generating  the  rudder  inputs  for 
producing  the  path  of  travel  desired. 

Since  the  course  of  the  sea  vessel  is  driven  by  the  rudder  inputs,  it  is 
necessary  to  compute  the  duration  of  the  rudder  inputs  as  the  sea  vessel  goes 
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through  each  turn.  This  can  be  easily  done  by  considering  the  angle  that  the  ship 
turns  through,  since  we  have  already  determined  the  time  it  takes  for  the  ship  to 
turn  through  a  complete  round.  Referring  to  the  angle  9  in  Figure  18,  the  angle 
that  the  sea  vessel  turns  through  can  be  computed  by  computing  180 °-9, 
based  on  the  geometrical  relationship  shown  in  Figure  38. 


Figure  38.  Geometrical  Relationship  for  Sea  Vessel  Turn  Angle 

The  turn  angle  computation  allows  us  to  compute  the  turn  duration  for 
every  turn  in  the  path.  Putting  all  the  turn  durations  into  a  combined  rudder  input 
for  the  simulation,  the  control  to  the  system  is  presented  in  Figure  39,  and  the 
resulting  output  corresponding  to  the  input  is  shown  in  Figure  40. 
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Rudder  Angle 


Figure  39.  Sea  Vessel  Rudder  Angle  Inputs 


Course  of  Sea  Vessel  on  X-Y  Plane 


Figure  40.  Sea  Vessel  Simulated  Course  of  Travel 
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Rudder  Angle 


From  the  results  obtained,  it  can  be  clearly  seen  that  the  sea  vessel 
responds  as  per  intended,  turning  with  the  changes  in  the  rudder  angles,  moving 
the  sea  vessel  along  the  sea  route  that  was  planned  and  discussed  in  the 
previous  chapter. 

Having  achieved  the  desired  path  of  travel  for  the  sea  vessel,  and  with  the 
associated  velocity  profile  generated  by  the  simulation  model,  we  are  now  ready 
to  convert  the  position  and  velocity  profile  into  the  TSPI  format  required  for  the 
weapon  IFS  model.  A  copy  of  the  Simulink  model  for  the  sea  route  simulation 
can  be  found  in  Appendix  C. 
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IV.  TIME  DISCRETIZATION 


A.  TRAVELING  TIME  COMPUTATION 

Having  computed  the  position  and  velocity  profiles  of  both  the  land  and 
sea  vehicle  at  every  point  on  their  paths  of  travel,  we  can  then  compute  the  time 
duration  for  the  vehicle  to  travel  between  each  consecutive  pair  of  points  on  the 
path.  This  can  be  done  using  the  kinematic  model  with  linear  acceleration,  as 
shown  in  Figure  41 . 


Figure  41 .  Kinematic  Model  for  Time  Duration  Computation 


With  the  distance  between  each  pair  of  consecutive  points  known,  the 
time  duration  can  be  computed  using  the  following  equation. 


At, 


->/+ 1 


uM+Ui 


(4.1) 
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After  we  have  computed  the  time  duration  the  vehicle  spent  traveling 
between  points,  we  can  then  proceed  to  recreate  the  path  of  travel  with  points  at 
discrete  time  intervals  to  generate  the  data  for  the  TSPI  inputs. 

B.  CONVERSION  TO  FIXED-TIME  INTERVAL 

To  recreate  the  path  of  travel  with  points  at  discrete  time  interval,  the 
approach  is  to  iteratively  step  through  each  point  in  the  path,  sum  up  the  time 
duration  for  the  vehicle  to  travel  to  the  current  point,  and  check  whether  the  sum 
exceeds  multiples  of  one  second.  If  the  total  time  duration  is  less  than  one 
second,  the  model  will  proceed  with  the  next  iterative  step  to  the  next  point  on 
the  path.  If  the  total  time  duration  is  more  than  one  second,  a  new  point  on  the 
new  path  with  a  discrete  time  interval  is  created  by  interpolating  between  the 
current  point  and  the  previous  point  for  which  the  time  duration  is  exactly  one 
second,  and  then  subtracting  one  second  from  the  total  time  duration.  This 
continues  until  the  total  time  duration  is  less  than  one  second,  and  the  model  will 
proceed  with  the  next  iterative  step  to  the  next  point  on  the  path.  A  visual 
representation  of  the  iterative  process  is  presented  in  Figure  42. 


Figure  42.  Creation  of  New  Points  at  Discrete  Time  Intervals 
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Besides  creating  the  new  points  at  every  discrete  time  interval,  it  is 
important  to  compute  the  velocity  of  the  vehicle  at  each  of  the  newly  created 
points  as  well.  This  can  be  done  during  the  same  iterative  process  by 
interpolating  the  velocity  of  the  new  point  from  the  velocities  of  the  vehicle  at  the 
points  between  which  the  new  point  is  being  created,  based  on  the  position  of  the 
new  point  relative  to  the  positions  of  the  points  for  the  interpolation.  This 
approach  is  acceptable  due  to  the  linear  acceleration  model  that  has  been 
adopted  for  the  model  earlier.  The  general  equation  for  the  interpolation  of  the 
velocity  of  the  new  point  j  from  those  of  point  i  and  (i  + 1)  is  given  in 
Equation  (4.2). 


Uj  =  ut  + 


~h 


(4.2) 


Having  derived  the  position  and  velocity  profile  of  the  new  points  for  the 
path  at  each  discrete  time  interval,  the  required  inputs  for  the  TSPI  file  can  now 
be  generated. 


C.  NEW  DATA  GENERATION  FOR  LAND  ROUTES 

Using  the  approach  that  has  been  discussed,  the  routes  that  we  have 
generated  previously  can  now  be  converted  to  ones  with  the  discrete  time 
intervals.  The  new  paths  for  the  routes  of  road  type  1  to  road  type  6  are  created 
and  presented  in  Figures  43  to  48,  respectively. 
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Figure  43.  Type  1  Route  at  Discrete  Time  Intervals 
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Figure  44.  Type  2  Route  at  Discrete  Time  Intervals 
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Figure  45.  Type  3  Route  at  Discrete  Time  Intervals 


Type  4  Road 


Figure  46.  Type  4  Route  at  Discrete  Time  Intervals 
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Figure  47.  Type  5  Route  at  Discrete  Time  Intervals 


Figure  48.  Type  6  Route  at  Discrete  Time  Intervals 
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Due  to  the  large  number  of  data  points  used  for  the  generation  of  the  six 
routes  for  type  1  to  type  6,  only  5  percent  of  the  points  are  presented  in  order  to 
see  the  new  data  points  on  the  plots  more  visibly.  Plotting  the  entire  range  of 
data  points  will  only  generate  plots  that  look  similar  to  those  produced  initially. 

From  the  plots  obtained,  it  can  be  observed  that  the  points  get  closer 
together  along  a  curvy  portion  of  the  route,  and  further  apart  when  the  route  gets 
straighter.  As  the  vehicle  travels  faster  along  a  straight  road,  the  amount  of 
distance  covered  with  each  second  of  time  is  greater.  Similarly,  when  the  road 
gets  curvier,  the  vehicle  speed  is  reduced,  and  the  amount  of  distance  that  can 
be  covered  in  a  second  of  time  is  reduced  accordingly.  As  such,  the  results 
obtained  are  as  expected;  therefore,  the  results  are  deemed  valid  and 
acceptable. 

A  copy  of  the  MATLAB  codes  for  the  computation  of  the  velocity  profile 
and  the  conversion  of  the  vehicle  position  along  the  route  to  discrete  time 
intervals  can  be  found  in  Appendix  D. 

D.  NEW  DATA  GENERATION  FOR  SEA  ROUTES 

As  the  number  of  points  for  the  path  of  the  sea  vessel  is  significantly 
lower,  the  plots  of  the  original  path  and  those  of  the  discrete  time  interval  are 
presented  in  Figures  49  and  50,  respectively,  for  ease  of  comparison. 
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Figure  49.  Original  Sea  Vessel  Route 


Course  of  Sea  Vessel  on  X-Y  Plane  -  Is  interval 


Figure  50.  Sea  Vessel  Route  at  Discrete  Time  Interval 
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Rudder  Angle  /  degrees  Rudder  Angle 


Similar  to  the  results  obtained  for  the  land  vehicle,  the  data  points  get 
closer  together  as  the  curvature  of  the  path  gets  higher,  and  further  apart  as  the 
curvature  of  the  path  decreases.  This  is  again  in  accordance  with  the 
expectations,  and  hence  the  results  are  valid  and  acceptable.  A  copy  of  the 
MATLAB  codes  for  the  conversion  of  the  sea  routes  to  discrete  time  intervals  can 
be  found  in  Appendix  E. 
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V.  ROUTE  DIGITIZATION 


The  past  few  chapters,  which  constitute  the  first  part  of  the  thesis,  focused 
on  generating  the  routes  of  interest  for  weapon  assessment  study  and  converting 
them  into  the  TSPI  inputs  that  can  be  used  in  the  weapon  IFS  model  for  analysis. 
For  the  rest  of  the  chapters,  which  constitute  the  second  part  of  the  thesis,  the 
aim  is  to  come  up  with  a  methodology  to  convert  an  actual  route  selected  from  a 
map  and  convert  it  into  a  digital  form  for  route  characterization  and  weapon 
assessment.  An  additional  requirement  for  this  methodology  is  a  Visual  Basic  for 
Applications  (VBA)  environment.  This  environment  enables  intended  users  to 
utilize  the  program  in  a  restricted  work  area  with  no  access  to  advanced  software 
for  performing  such  analysis. 

A.  SELECTING  THE  ROUTE  OF  INTEREST 

For  the  route  selection,  we  want  to  take  an  approach  that  is  simple,  user- 
friendly,  and  easily  accessible.  With  that  in  mind,  we  start  off  with  the  selection  of 
the  map  route  using  Google  Maps  [15],  a  web-based  electronic  map,  commonly 
used  by  many  around  the  world.  Google  Maps  is  easily  accessed  using  web 
browsers  and  an  active  internet  connection.  The  limitation  to  this  approach  is  that 
the  routes  available  for  selection  are  limited  to  areas  that  Google  has  managed 
to  access  and  map  out. 

Selection  of  a  route  using  Google  Maps  [15]  is  simple  and  straightforward, 
as  this  can  be  done  using  the  “Directions”  function  that  can  be  accessed  with  just 
a  click  of  a  mouse  button.  Thereafter,  the  selection  of  the  route  of  interest  can  be 
done  by  clicking  on  the  locations  on  the  map  to  specify  the  start  and  end  points 
of  the  route.  Google  Maps  [15]  will  then  process  and  generate  the  possible 
routes  to  get  from  the  start  point  to  the  end  point.  Then,  users  can  choose  from 
the  routes  that  Google  Maps  [15]  has  suggested  or  can  specify  the  route  that 
they  want  by  including  additional  waypoints.  An  illustration  of  this  for  a  route 
along  the  Carmel  Valley  road  is  presented  in  Figure  51 . 
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Figure  51 .  Selection  of  Route  of  Interest 


After  selecting  the  route  of  interest,  we  can  then  proceed  to  extract  the 
route  information  and  convert  it  into  a  digital  form  that  allows  for  further 
processing  of  the  information  that  we  are  keen  to  sieve  out. 

B.  IMPORTING  THE  ROUTE  OF  INTEREST 

One  of  the  common  file  formats  for  presenting  map  information  is  the  GPX 
file  format  that  GPS  devices  utilize.  For  the  purpose  of  this  thesis,  we  use  the 
GPX  file  format  for  the  extraction  of  the  route  information.  While  Google 
Maps  [15]  provides  a  user-friendly  interface  for  route  selection,  it  does  not  allow 
for  the  extraction  of  the  route  information.  As  such,  we  need  to  utilize  other 
avenues  to  get  the  data  we  need  from  Google  Maps.  A  quick  search  on  the  web 
will  reveal  many  tools,  available  for  purchase  or  at  no  cost,  for  performing  this 
task.  With  the  key  consideration  of  accessibility  in  mind,  a  web-based  conversion 
tool  is  preferred.  Therefore,  for  this  study,  a  third-party  website,  MapsToGPX 
[16],  is  chosen  not  only  for  its  accessibility  online  but  its  simplicity  of  use. 
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Extraction  of  the  map  information  can  be  easily  done  by  inserting  the 
Uniform  Resource  Locator  (URL),  also  commonly  known  as  the  web  address,  of 
the  route  that  we  have  selected  in  Google  Maps  [15]  earlier  and  clicking  the 
“Let’s  Go”  button  to  process  the  route  information  into  a  downloadable  GPX  file. 
An  illustration  of  the  usage  of  the  map  conversion  tool  is  provided  in  Figure  52. 
With  this  file,  the  information  of  all  the  points  that  define  the  route  selected  in 
Google  Maps  is  retrieved  in  digital  format. 


Figure  52.  Generation  of  Route  Data  in  GPX  Format 


After  obtaining  the  route  information  in  GPX  file  format,  we  must  import 
the  information  into  a  software  tool  so  that  the  information  can  be  read  and 
processed  easily.  For  this  purpose,  both  MATLAB  and  Microsoft  EXCEL  are 
capable  of  importing  the  GPX  file  and  presenting  the  map  data  in  terms  of  their 
latitude  and  longitude  coordinates  in  array  form.  With  this,  we  are  ready  to 
process  the  information  for  further  analysis  of  the  route. 
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VI.  ROUTE  PROCESSING  AND  ANALYSIS 


Using  the  approach  discussed  earlier  under  “Map  Coordinates 
Transformation”  in  Chapter  I,  we  can  convert  the  map  data  extracted  from  the 
GPX  file  from  the  geodetic  latitude  and  longitude  coordinates  to  the  ENU 
coordinate  system.  In  this  way,  the  Easting  and  Northing  coordinates  are 
equivalent  to  the  familiar  Cartesian  coordinate  system. 

The  key  parameter  that  we  want  to  extract  from  the  route  is  the  standard 
deviation  of  the  change  in  heading  of  the  route.  This  parameter  enables  us  to 
classify  the  route  as  one  of  the  six  road  types  presented  in  Table  2  in  Chapter  II. 
To  compute  the  standard  deviation  of  the  change  in  heading,  first  we  need  to 
standardize  the  route  make-up  for  a  fair  statistical  computation  of  the  heading 
and  change  in  heading.  This  is  done  by  converting  the  route  from  one  comprising 
irregular  route  points  to  one  with  regularly  spaced  route  points,  and  for  this 
thesis,  the  route  points  are  spaced  ten  meters  apart. 

A.  CONVERSION  TO  REGULAR  INTERVALS 

For  the  conversion  of  the  route  points  to  regular  intervals,  the  approach  is 
similar  to  that  used  for  conversion  to  the  fixed-time  interval  that  we  discussed  in 
Chapter  IV. 

To  recreate  the  path  of  travel  with  points  at  equal  intervals  of  ten  meters 
apart,  the  process  required  is  to  iteratively  step  through  each  point  in  the  path, 
sum  up  the  distance  traveled  for  the  vehicle  to  reach  the  current  point,  and  check 
whether  the  sum  exceeds  multiples  of  ten  meters.  If  the  total  distance  traveled  is 
less  than  ten  meters,  the  model  will  proceed  with  the  next  iterative  step  to  the 
next  point  on  the  path.  If  the  total  time  duration  is  more  than  ten  meters,  a  new 
point  on  the  new  path  with  regularly  spaced  points  is  created  by  interpolating  for 
the  position  between  the  current  point  and  the  previous  point  for  which  the 
distance  traveled  is  exactly  ten  meters,  and  then  subtracting  ten  meters  from  the 
total  distance  traveled.  This  continues  until  the  total  distance  travelled  is  less 
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than  ten  meters,  and  the  model  will  proceed  with  the  next  iterative  step  to  the 
next  point  on  the  path.  A  visual  representation  of  the  iterative  process  is 
presented  in  Figure  53. 


(x5,y5) 


Figure  53.  Creation  of  New  Points  at  Regular  Distance  Interval 


The  general  equation  for  the  interpolation  of  the  x  and  y  coordinates  of  the 
new  point  j  from  those  of  point  i  and  (i  +  1)  is  given  as 


AS; 


Xj=Xi+- 


Ujixi+ i~xi) 


si+l  ~si 


(6.1) 


yj  =  T;  + 


A^.(y,+1-y,) 


SM~St 


(6.2) 


This  approach  is  used  to  convert  the  Carmel  Valley  route  extracted  from 
Google  Maps  [15]  into  a  route  made  up  of  points  spaced  ten  meters  apart.  The 
original  and  newly  converted  routes  are  presented  in  Figures  53  and  54, 
respectively. 
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Figure  54.  Original  Route  Points  Imported 


Figure  55.  New  Regularly  Spaced  Route  Points 
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From  the  comparison  of  the  visual  plots  of  the  two  routes,  we  can  see  that 
the  original  route  has  been  successfully  converted  into  one  with  the  route  points 
regularly  spaced  at  ten  meters  apart.  Having  established  the  methodology  for 
standardizing  the  route  in  terms  of  their  make-up,  we  can  then  proceed  to 
conduct  the  statistical  analysis  for  the  converted  route. 

B.  CHANGE  OF  HEADING  COMPUTATION 

To  compute  the  heading,  we  need  to  establish  the  heading  of  each  of  the 
points  on  the  path.  We  do  this  by  calculating  the  slope  of  the  line  joining  the 
points  immediately  before  and  after  the  point  of  interest.  After  the  standardization 
of  the  route  points  earlier,  the  regularly-spaced  route  points  allow  this  approach 
to  give  a  good  approximation  for  the  slope  of  the  tangent  to  the  route  at  the 
current  point.  This  methodology  is  applied  for  the  computation  of  the  slope  at 
every  point  along  the  route.  The  computation  of  the  slope  of  the  tangent  to  the 
route  and  the  heading  at  a  route  point  is  presented  in  Figure  56. 
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Figure  56.  Computation  of  Heading 

With  this  approach,  the  heading  of  every  point,  less  the  start  and  end 
points,  along  the  route  can  be  computed.  With  the  heading  on  each  point 
established,  the  change  of  heading  can  then  be  easily  computed  by  taking  the 
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difference  between  the  heading  of  the  current  point  and  that  of  the  point  prior. 
The  computation  of  the  change  in  heading  of  the  point  is  illustrated  in  Figure  57. 
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Figure  57.  Computation  of  Change  in  Heading 

The  approach  is  tested  on  a  short  winding  route,  shown  in  Figure  58.  The 
results  of  the  heading  computations  are  presented  in  Figures  59  and  60. 
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Figure  59.  Heading  of  Points  on  Test  Route 
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Figure  60.  Change  of  Heading  of  Points  on  Test  Route 
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With  the  positive  results  seen  from  the  test,  we  are  able  to  proceed  to 
compute  the  change  in  heading  for  all  the  points  on  the  actual  route.  The  newly 
computed  data  for  change  in  heading  then  allows  us  to  carry  out  statistical 
analysis  to  calculate  the  standard  deviation  of  the  distribution.  Using  the  same 
example  of  the  Carmel  Valley  road,  we  find  the  standard  deviation  of  the  change 
in  heading  was  calculated  to  be  5.643  degrees,  which  classifies  the  route  as  a 
type  3  road.  This  is  close  to  the  results  obtained  by  Driels  [10],  as  discussed  in 
Chapter  II,  and  helps  to  validate  the  approach  adopted  in  this  study.  Having 
validated  the  approach  by  obtaining  similar  results  with  previous  works,  we  then 
proceed  to  apply  the  program  to  routes  in  other  countries  where  mapping  of 
routes  may  be  less  extensive  than  that  in  the  United  States.  The  program  is  used 
for  selected  routes  in  North  Korea,  Syria,  and  Afghanistan.  The  program  works 
well,  and  the  results  are  presented  in  Figures  61  to  Figure  63,  respectively. 


Figure  61 .  Processed  North  Korean  Route — Type  2 
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Figure  62.  Processed  Syrian  Route — Type  2 


Figure  63.  Processed  Afghanistan  Route — Type  3 
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From  the  results  for  the  routes  shown,  it  can  be  seen  that  the  approach 
works  for  the  routes  that  can  be  selected  using  Google  Maps  and  imported  into 
the  software  tool. 

The  model  is  initially  developed  and  tested  in  MATLAB.  As  discussed 
earlier,  a  fresh  copy  is  translated  into  the  VBA  environment  to  suit  the  intended 
user.  With  the  click  of  a  single  macro  button  by  the  user,  the  VBA  program  will 
automatically  extract  the  route  data  from  the  GPX  file  and  analyze  the  route  to 
compute  the  standard  deviation  in  the  change  of  heading,  thereby  categorizing 
the  route  into  the  appropriate  road  type  based  on  the  results  obtained.  The  user 
interface  in  the  VBA  environment  is  presented  in  Figure  64. 
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Download  the  gpx  file  to  the  following  folder  C:\Weapon\GPXfile  and  save  the  file  as  "Route.gpx" 
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Click  this  button  - 
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Click  me! 
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Step  6: 

Go  back  to  Step  1  for  new  route,  or  exit  this  program  if  you  are  done. 
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Figure  64.  Route  Analysis  in  VBA  Environment 


Both  models  are  able  to  produce  the  same  results,  verifying  that  the 
working  approach  is  successfully  implemented  in  an  interface  suitable  for  the 
intended  users. 

A  copy  of  the  code  for  the  route  analysis  implemented  for  MATLAB  and 
VBA  can  be  found  in  Appendix  F  and  Appendix  G,  respectively. 


71 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


72 


VII.  CONCLUSION  AND  RECOMMENDATIONS 


A.  CONCLUSION 

This  thesis  has  achieved  the  objectives  set  by  successfully  generating  the 
TSPI  files  for  different  synthesized  routes  that  are  representative  of  actual  target 
routes  of  interest.  With  these  TSPI  files  obtained,  the  routes  can  then  be  fed  into 
the  weapon  IFS  model  for  weapon  assessment  against  a  moving  target  traveling 
along  the  specific  route  type.  The  results  of  the  weapon  assessment  can  then  be 
compiled  to  create  a  quick  look-up  reference  for  assessment  of  the  effectiveness 
of  a  specific  weapon  against  moving  targets  moving  along  route  of  a  particular 
road  type. 

A  program,  written  in  the  VBA  environment,  is  also  developed  to  enable 
quick  categorization  of  an  actual  road  of  interest  into  one  of  the  road  profiles 
described,  and  this  program  enables  the  weaponeer  to  obtain  quick  estimates  of 
the  effectiveness  of  various  weapon  types  available  and  facilitates  a  quick 
selection  of  weapon  type  for  optimal  probability  of  kill. 

The  key  products  of  this  thesis  include: 

•  the  MATLAB  program  for  generation  of  synthesized  constrained 
routes; 

•  the  MATLAB  program  for  generating  TSPI  files  for  synthesized 
routes.; 

•  the  Simulink  program  for  generating  expected  unconstrained 
routes  for  sea  vessels; 

•  the  MATLAB  program  for  generating  TSPI  files  for  sea  routes; 

•  the  MATLAB  program  for  digitization  of  map  routes  and 
categorization  of  map  routes  into  specific  road  types;  and 
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•  the  VBA  program  for  digitization  of  map  routes  and  categorization 
of  map  routes  into  specific  road  types. 

B.  RECOMMENDATIONS 

This  thesis  can  be  further  extended  into  the  following  areas  for  future 

works: 

•  The  velocity  computation  can  be  done  using  the  Simulink  model 
discussed  earlier  to  better  capture  the  effects  of  vehicle  dynamics 
to  produce  results  that  are  more  representative  of  the  actual  vehicle 
performance. 

•  The  velocity  computation  can  be  optimized  to  determine  the  optimal 
time  for  a  given  vehicle  to  travel  a  route  in  minimal  time. 

•  The  TSPI  data  obtained  can  be  compared  to  actual  vehicle  trials  to 
compare  the  theoretical  and  empirical  results  to  verify  the  approach 
and  assumptions  used. 
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APPENDIX  A.  MATLAB  CODE  FOR  SYNTHESIZED  ROUTE 

GENERATION 


clear  all 
clc 

heading_limit=2 ;  %  #  of  sigmas  heading  is 

limited 

% - New  data  set  for  road  types  1-6 - 

mean_heading=  [135  135  135  135  135  135]; 
sigma_heading=  [25  29  41  53  66  73]; 
mean_dheading=  [0  0  0  0  0  0]; 
sigma_dheading=  [0.5  2.5  7  12  17  25]; 
total_km  =  [40  40  40  40  40  40]; 

road_type=menu ( '  Choose  a  road  type  'Type  1', '  Type 

2  ',.... 

'  Type  3  ' , '  Type  4  ' , ' Type  5 ' , '  Type  6  ' ) ; 

r_type= [ ' Type  1 ' ; ' Type  2 ' ; ' Type  3 ' ; ' Type  4 ' ; ' Type 
5 ' ; ' Type  6 ' ] ; 

r_type_data=cellstr (r_type) ; 

%  get  statistics  for  the  chosen  road  type: 

av_heading  =  90  -  mean_heading (road_type) ;  %  rotate 

coordinates  from  north  to  east 

sig_heading=sigma_heading (road_type) ; 

dhmean  =  mean_dheading (road_type) ; 

dhstd_dev  =  sigma_dheading (road_type) ; 

h_limit  =  heading_limit *sig_heading;  %  heading 

constraint  limit 


heading (1)  =  av_heading;  %  starting  heading 

seg_dist  =10;  %  10  meter  segments 

npts  =  tot al_km ( road_type )  *  1000  /  seg_dist; 


%  Randomly  sample  from  the  normal  distribution 
%  Each  r  value  gives  you  the  change  in  heading  for  the 
next  10  m  of  road 
for  i=l : npts 

%r(i)  =  random (' norm dhmean, dhstd_dev) ; 
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r (i) =randn*dhstd_dev+dhmean; 
end 

stdev_data=std (r )  %  Calculate  the  actual  sigma  of 

this  data  set 


o 

o 

%  Here's  the  road  building  part,  which  turns  a  heading 
change  into  x  and  y 

%  coordinates  starting  at  easting(l),  northing(l). 
easting(l)  =  615788; 
northing(l)  =  4035674; 

for  i=2 : npt s+1 


heading (1)  =  heading (l-l)  +  r(i-l); 

%try  to  keep  road  from  making  circles: 


if  heading  (1)  >  (av_headmg 

(av_heading  -  h_limit) 

heading (i)  =  heading (i-1) 

end 

easting(i)  =  easting(i-l)  + 
1) *pi/180) ; 

northing (i)  =  northing ( i-1 ) 

1) *pi/180) ; 
end 


+  h_limit)  |  heading (i)  < 

-  r  (i-1) ; 

seg_dist*cos (heading (i- 
+  seg_dist*sin (heading (i- 


%  plot  the  road 
figure (2 ) 

plot (easting/1000 ,  northing/1000, ' LineWidth ' , 1 . 25 ) 

axis ( ' equal ' ) ;  grid 

%  title (' Synthesized  Road') 

title (r_type_data (road_type) ) 

xlabel ( 'Easting  (km) ' ) 

ylabel ( 'Northing  (km) ' ) 
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APPENDIX  B.  SIMULINK  MODEL  FOR  VELOCITY  COMPUTATION 
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APPENDIX  C.  SIMULINK  MODEL  FOR  SEA  ROUTE  GENERATION 


Subsystem 
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APPENDIX  D.  MATLAB  CODE  FOR  SYNTHESIZED  ROUTE 
PROCESSING  INTO  TSPI  FORMAT 


clear  all 
clc 

load ( ' Type3 . mat ' ) 
for  i=l : length (type3 ) 

newroute (i, 1)  =  type3(i,l)  -  type3(l,l);  %  initialise 
1st  data  point  as  origin 

newroute ( i , 2 )  =  type3(i,2)  -  type3(l,2); 
end 

%%  Initialization  for  Computing  Speed 
veh_limit  =60;  %  vehicle  speed  limit  in  mph 
m2  ft  =  3.28084;  %  conversion  from  metres  to  feet 
mph2ms  =  1.60934  *  1000  /  3600;  %  conversion  from  mph 
to  m/s 

coef f_f rict ion  =  0.7;  %  coefficient  of  friction 
safety_f actor  =  1;  %  speed  safety  factor 
interval  =  10;  %  distance  between  subsequent  points  is 

1 0m 

check  =10;  %  max  number  of  points  to  check  prior  and 

after 

decel_limit  =0.6  *  9.81;  %  deceleration  limit,  from  g 
to  m/s2 

accel_limit  =0.3  *  9.81;  %  acceleration  limit,  from  g 
to  m/s2 

newroute ( 1 , 3 )  =  0;  %  initialize  distance 
newroute (1,4:9)  =  0 ; 

newroute ( 1 , 1 0 )  =  safety_f actor  *  min(sqrt(15  *  lel5  * 
m2ft)  *  (coef f_f riction  +  0.01  *  0 ) , veh_limit )  * 

mph2ms;  %  initialize  vehicle  at  speed  limit,  with  bank 
angle  =  0 

for  i=2 : length (newroute ( : , 1) ) -1 

newroute ( i , 3 )  =  ( (newroute (i, 1 ) -newroute (i-1,  1 )). A2  + 

(newroute (i, 2) -newroute (i-1, 2) ). A2)  .A  0.5;  %  there  are 
values  less  than  zero  as  previous  computation  was  based 
on  accumulated  distance  travelled 
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newroute (i, 4 )  =  (newroute ( i , 2 )  -  newroute (i-1 , 2 ) )  / 

(newroute ( i, 1 )  -  newroute ( i-1 , 1 )) ;  %  gradient  between 

point  i  and  i-1 

newroute (i, 5)  =  (newroute ( i+1 , 2 )  -  newroute ( i , 2 ) )  / 

(newroute ( i+1 , 1 )  -  newroute (i, 1 )) ;  %  gradient  between 
point  i  and  i+1 

newroute (i, 6)  =  (newroute ( i , 4 )  *  newroute ( i , 5 )  * 

(newroute ( i+1 , 2 )  -  newroute ( i-1 , 2 ) )  +  newroute ( i , 4 )  * 

(newroute ( i, 1 )  +  newroute (i+1, 1) )  -  newroute (i, 5)  * 

(newroute ( i, 1 )  +  newroute ( i-1 , 1 )) )  /  (2  * 

(newroute ( i, 4 )  -  newroute (i, 5) )) ;  %  compute  x- 

coordinate  of  centre  of  circle  of  curvature 

newroute (i, 7 )  =  -  1  /  newroute ( i , 5 )  *  (newroute ( i , 6 ) 

-  (newroute (i, 1 )  +  newroute ( i+1 , 1 ) )  /  2)  + 

(newroute ( i, 2 )  +  newroute ( i+1 , 2 ) )  /  2;  %  compute  y- 
coordinate  of  centre  of  circle  of  curvature 

newroute (i, 8)  =  min (( (newroute ( i , 6 )  -  newroute ( i , 1 ) ) 

A  2  +  (newroute (i, 7 )  -  newroute ( i, 2 ) )  A  2)  A  0.5,lel5); 

%  radius  of  curvature  in  metres,  lel5  for  those  with 
infinite  values,  i.e  collinear  points 
newroute ( i , 9 )  =0;  %  road  bank  angle 
newroute (i, 10 )  =  safety_f actor  *  min(sqrt(15  * 
(newroute ( i, 8 )  *  m2ft)  *  (coef f_f riction  +  0.01  * 

newroute ( i , 9 ))), veh_limit )  *  mph2ms;  %  max  speed  limit 

at  point  i  in  m/s,  converted  from  mph  using  radius  of 
curvature  in  feet 

newroute ( i , 1 6 )  =  sum (newroute ( 1 : i, 3 )) ;  %  Displacement 
end 

%%  Generate  Speeds  at  points 

newroute ( 1 , 1 1 )  =  sqrt (newroute (2, 10) A2  +  2  * 
decel_limit  *  interval) ; 

newroute ( 1 , 12 )  =  sqrt (newroute (2 , 1 0 ) A2  -  2  * 

accel_limit  *  interval) ; 

newroute ( 1 , 13 )  =  max (newroute ( 1 , 11 ), newroute ( 1 , 12 )) ; 

for  i=2 : length (newroute ( : , 1) ) -1 

%  check  for  deceleration  constraint 
acc_checkl  =  NaN (check, 1 ) ; 
vel_checkl  =  NaN (check, 1 ) ; 

for  j=l :min ( (length (newroute) -1)  -  i,  check)  % 

iteration  for  check  on  speed  limits  on  points  ahead 
if  newroute (i, 10 )  >  newroute (i+j , 10 ) 
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vel_checkl ( j )  =  sqrt (newroute (i+ j , 10 ) A2  +  2  *  j  * 
decel_limit  *  interval) ; 
else 
end 
end 

newroute (i, 11 )  =  min (min (vel_checkl ), newroute (i, 10 )) ; 

%  check  for  acceleration  constraint 
acc_check2  =  NaN (check, 1 ) ; 
vel_check2  =  NaN (check, 1 ) ; 

for  j=l:min(i-l,  check)  %  iteration  for  check  on 
speed  limits  on  points  prior 

if  newroute (i, 10 )  >  newroute (i-j , 10 ) 

vel_check2 ( j )  =  sqrt (newroute (i-j , 10 ) A2  +  2  *  j  * 

accel_limit  *  interval) ; 
else 
end 
end 

newroute ( i , 12 )  =  min (min (vel_check2 ), newroute (i, 10 )) ; 
newroute ( i , 13 )  =  min (newroute ( i, 11 ), newroute ( i , 12 )) ; 

%  Actual  speed/  highest  speed  possible 
end 

%%  Generate  time  between  points 
for  i=l : length (newroute (:, 1 )) -2 

newroute (i, 14 )  =  (newroute (i+1, 13) A2  - 
newroute ( i , 13 ) A2 )  /  (2  *  interval);  %  Acceleration 

newroute ( i , 15 )  =  2  *  interval  /  (newroute (i, 13)  + 
newroute (i+1, 13) ) ;  %  time  from  point  i  and  i+1 
end 

%%  Generate  Points  at  1  sec  Intervals 
carryover  =  0; 
k  =  2 ; 

time  =  newroute (:, 15 ) ; 
newroutetime ( 1 , 1 )  =  newroute (1, 1) ; 

newroutetime ( 1 , 2 )  =  newroute ( 1 , 2 ) ; 

newroutetime ( 1 , 4 )  =  0; 

for  j=l : length (time) -1 

if  time(j)  +  carryover  <  1 

carryover  =  time  (j)  +  carryover; 

else 
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for  count=l : floor ( (carryover  +  time(j))/l) 
newroutetime (k, 8 )  =  k  -  1;  %  total  time 
newroutetime (k, 1 ) =newroute (  j  ,  1 )  + (newroute ( j+1 , 1 ) - 
newroute ( j , 1 ) ) * (count-carryover) /time ( j ) ; 

newroute time (k, 2 ) =newroute (  j  ,  2  )  + (newroute (j+1, 2) — 
newroute ( j , 2 ) ) * (count-carryover) /time ( j ) ; 

newroutetime (k, 9) =newroute ( j , 13 ) + (newroute (j+1, 13 ) - 
newroute ( j , 13) )* (count-carryover ) /time ( j ); %  speed  at 
point  k 

k=k+l ; 
end 

carryover  =  carryover  +  time(j)  -  count; 
end 
end 

%%  Heading  and  Speed  for  Time-Based  route 
for  i=2 : length (newroutetime) -1  %  heading 
newroutetime (i, 4 )  =  newroutetime ( i+1 , 1 )  - 
newroutetime ( i-1 , 1 ) ;  %  delta  x 

newroutetime ( i , 5 )  =  newroutetime ( i+1 , 2 )  - 

newroutetime ( i-1 ,  2 )  ;  %  delta  y 

newroutetime (i, 6)  = 

atan2 (newroutetime (i, 5) , newroutetime (i, 4 )) ;  %  heading 

with  horizontal  as  0  rad 

if  newroutetime (i, 6)  >  pi/2  &&  newroutedist (i, 6)  <  pi 

%  heading  with  vertical  as  0  rad 

newroutetime ( i , 7 )  =  2  *  pi  +  (  pi  /  2  - 

newroutetime (i, 6) ) ; 
else 

newroutetime ( i , 7 )  =  pi  /  2  -  newroutetime ( i , 6 ) ; 
end 

newroutetime (i, 10)  =  newroutetime (i, 9)  * 

sin (newroutetime (i, 7 ) ) ; 

newroutetime (i, 11)  =  newroutetime (i, 9)  * 

cos (newroutetime (i, 7) ) ; 
end 

%%  Conversion  From  ENU  to  ECEF  to  Geodetic 

a  =  6378137.0;  %  WGS-84  Earth  semimajor  axis  (m) 

b  =  6356752.3142;  %  WGS-84  Earth  semiminor  axis  (m) 

f  =  (a-b)/a;  %  Ellipsoid  Flatness 

e_sq  =  f  *  (2  -  f ) ;  %  Square  of  Eccentricity 
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epsilon  =  e_sq  /  (1  -  e_sq) ; 

h  =  0; 


routel  (1,1)  = 
point  lat/lon 
routel  (1,2)  = 
routel (1,3) 
routel (1,4) 
routel (1,3) 
axis  along 
routel (1,5) 


0; 

to 

0; 


%  initialise  synthesized  route  start 

0,  0 


=  sin (routel (1, 1) ) ; 


sin  (lat) 

=  a  /  sqrt(l-e_sq  *  routel (1,3)  * 

) ;  %  N,  distance  from  the  surface  to 
the  ellipsoid  normal 
=  (h  +  routel (1,4))  * 

cos  (routel ( 1 , 2 )) ;  %  converting  to 

routel (1,6)  =  (h  +  routel (1,4))  * 

sin  (routel (1,2) ) ; 

(h  + 


the  z- 


routel (1,7) 


cos (routel (1, 1) )  * 

ECEF  x-coordinates 
cos (routel  (1, 1) )  * 

converting  to  ECEF  y-coordinates 
(1  -  e_sq)  *  routel  (1,4))  * 


Is 


sin (routel ( 1 , 1 )) ;  %  converting  to  ECEF  z-coordinates 

for  i=l : length (newroutetime)  %  converting  points  for 
intervals  back  to  Lat/Lon 

route3(i,18)  =  -sin ( routel ( 1 , 2 ) )  *  newroutetime ( i , 1 ) 

-  sin (routel (1, 1) )  *  cos ( routel ( 1 , 2 ) )  * 

newroutetime ( i , 2 )  +  cos ( routel  ( 1 , 1 ) )  *  cos  ( routel  ( 1 , 2 ) ) 

*  newroutetime (i, 3)  +  routel (1,5);  %  convert  back  to 
ECEF  x-coordinate 

route3(i,19)  =  cos ( routel  ( 1 , 2 ) )  * 

sin  (routel ( 1 , 1 ) )  *  sin  (routel ( 1 , 2 ) ) 

*  sin (routel (1,2 
+  routel (1,6);%  c 


+  cos (routel (1,1 
newroutetime (i, 3 
y-coordinate 

route3  (i, 20 )  = 
sin (routel (1,1) ) 
convert  back  to 
route3  ( i , 2 1 )  = 

%  p 

route3 (i, 22 ) 

*  b)  )  ;  %  q 

route3  (i, 23 )  = 
(sin (route 3 (i, 22 
(cos (route3 (i, 22 
route3  (i, 24 )  = 
Lon 
end 


newroutetime ( i , 1 )  - 

*  newroutetime ( i , 2 ) 

\  * 


cos (routel  (1, 1) )  * 

*  newroutetime ( i, 3 
ECEF  z-coordinate 
=  sqrt ( route3 ( i , 1 8 ) A2 

=  atan ( (route3 (i, 20 ) 

=  atan2 ( ( route3 ( i , 2 0 ) 
) ) A3)  ,  (route3 (i, 2 

) ) A3) ) ;  %  Lat 
atan2  (route3 (i, 19) 


invert  back  t 

o 

ECEF 

newroutet ime 

(i 

,  2 )  + 

+  routel ( 1 , 

7) 

;  % 

:  +  route3 ( i , 

19 

)  A2)  ; 

*  a) /  (route 

3  ( 

i/  21) 

+  epsilon  * 

b 

* 

:1)  -  e_sq  * 

a 

* 

, route3 (i, 18 

)  ) 

•  £- 

r  ° 
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NRTLatLon ( : , 1) =route3 ( : , 23) *180/pi;  %converting  Is  x- 
coordinates  back  to  Lon 

NRTLatLon ( : , 2 ) =route3 ( : , 2 4 ) * 1 80 /pi ;  %converting  Is  y- 
coordinates  back  to  Lat 
%%  Extracting  TSPI 

TSP I  =  [newroutetime (:, 8)  NRTLatLon (:, 1 )  NRTLatLon (:, 2 ) 
newroutetime ( : , 10 )  newroutetime (  : , 11 ) ]  ; 
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APPENDIX  E.  MATLAB  CODE  FOR  SEA  ROUTE  PROCESSING 


clear  all 
clc 

%%  Simulate  the  Simulink  Model 

sim ( ' Sea_Vessel_Plant . six ' , 90 ) ;  %  runnning  Model  in 
Simulink 

%%  Read  the  timeseries  object  from  the  data  file 
load  boat . mat 
%%  Parcing 

t ime_sim=boat_result s . time ; 
route ( : , 1 ) =boat_result s . Data ( : , 1 ) ; 
route ( : , 2 ) =boat_result s . Data ( : , 2 ) ; 
rudder  =  boat_result s . Data ( : , 3 ) ; 
speed  =  boat_result s . Data ( : , 4 ) ; 

%%  Convert  Intervals  of  10m 

for  i=l : length (route (:, 1 )) -1  %finding  difference  in  x 

and  y  for  points  i  and  i+1 

dif f lx ( i ) =route ( i+1 , 1 ) -route ( i, 1 ) ; 

diff ly (i) =route (i+1, 2 ) -route (i, 2 ) ; 

end 

k=2 ;  %initiatise  k=2 

carryover (1) =0;  %  initiatise  carryover=0 
addx=0;  addy=0; 

newroute= [ route (1,1)  route (1,2) ] ; 
for  j  =  1 : i 

distance (j) = (difflx (j) A2+diffly (j) A2 ) A0 . 5;  % computing 
distance  between  points  j  and  j+1 
if  carryover (j)  >  0 

if  (carryover (j) +di stance (j) ) >=10 

for  count =1 : floor ( (carryover (j) +di stance (j) ) /10) 
if  count==l 

newroute (k, 1) =newroute (k- 

1. 1)  idifflx (j) /distance (j) * ( 1 0 -carryover ( j ) ) faddx; 

newroute (k, 2) =newroute (k- 

1.2)  +dif f ly ( j ) /distance ( j ) * (10-carryover ( j ) ) faddy; 

k=k+l; 

else 
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newroute (k, 1) =newroute (k- 

1.1)  +difflx ( j) /distance ( j) *10; 

newroute (k, 2) =newroute (k- 

1.2)  +diffly ( j) /distance ( j) *10; 

k=k+l; 

end 

end 

carryover ( j+1 ) =di stance ( j ) - (10-carryover ( j ) ) - 
(count-1) *10; 

addx=dif f lx ( j ) / distance ( j ) * carryover ( j+1 ) ; 
addy=dif f ly ( j ) / distance ( j ) * carryover ( j+1 ) ; 
else 

carryover ( j+1 ) =carryover ( j ) +di stance ( j ) ; 
addx=addx+dif f lx ( j ) ; 
addy=addy+dif f ly ( j ) ; 
end 
else 

if  distance ( j ) >=1 0  %  check  if  distance  to  next 
point  is  >  10m 

for  count =1 : floor (distance ( j ) /10) 
newroute (k, 1) =newroute (k- 

1.1)  +difflx ( j) /distance ( j) *10; 

newroute (k, 2) =newroute (k- 

1.2)  +diffly ( j) /distance ( j) *10; 

k=k+l; 

end 

carryover ( j+1 ) =di stance ( j ) - (count ) *1 0 ; 
addx=dif f lx ( j ) / distance ( j ) * carryover ( j+1 ) ; 
addy=dif f ly ( j ) / distance ( j ) * carryover ( j+1 ) ; 
elseif  distance ( j ) <10  %  check  if  distance  to  next 
point  is  <  10m 

carryover ( j  +  1 ) =di stance ( j )  ; 
addx=dif f lx ( j  )  ; 
addy=dif f ly ( j  )  ; 
end 
end 
end 

%%  Speed 

for  i=l : length (newroute (:,  1) ) 

row_index  =  find (newroute (i, 1 ) <route (:, 1 ), 1 ) -1 ; 
if  row_index  <  length (newroute (:, 1 ) ) 
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newroute (i, 3)  =  (newroute ( i , 1 )  -  route ( row_index, 1 ) ) 

/  (route (row_index+l , 1 )  -  route ( row_index, 1 ) )  * 

( speed (row_index  +  1)  -  speed (row_index) )  + 

speed ( row_index) ;  %  interpolate  speed 

end 
end 

%%  Generate  time  between  points 
interval  =  10; 

for  i  =  l : length (newroute (:,  1 )) -1 

newroute ( i , 4 )  =  2  *  interval  /  (newroute ( i, 3 )  + 

newroute ( i+1 , 3 )) ;  %  time  from  point  i  and  i+1 
end 

%%  Generate  Points  at  1  sec  Intervals 
carryover  =  0; 
k  =  2; 

time  =  newroute ( : , 4 ) ; 
newroutetime ( 1 , 1 )  =  time_sim(l); 

newroutetime ( 1 , 2 )  =  newroute (1, 1) ; 

newroutetime (1, 3)  =  newroute ( 1 , 2 ) ; 
for  j=l : length (time) -1 

if  time(j)  +  carryover  <  1 

carryover  =  time  (j)  +  carryover; 

else 

for  count=l : floor ( (carryover  +  time(j))/l) 
newroutetime (k, 1 )  =  k  -  1; 

newroutetime (k, 2 ) =newroute ( j , 1 )  +  (newroute ( j  +  1 ,  1 ) 
newroute ( j , 1 ) ) * (count-carryover) /time ( j ) ; 

newroute time (k, 3 ) =newroute ( j , 2 ) + (newroute ( j+1 , 2 ) 
newroute ( j , 2 ) ) * (count-carryover) /time ( j ) ; 
k=k+l; 
end 

carryover  =  carryover  +  time(j)  -  count; 
end 
end 

%%  Speed 

for  i  =  l:length (newroutetime ( :  ,  1 ) ) 

row_index  =  find (newroutetime (i, 2 ) <route (:, 1) , 1) -1; 
if  row_index  <  length ( route (:, 1 ) ) 
newroutetime ( i ,  4 )  =  (newroutetime (i, 2 )  - 

route ( row_index, 1 ) )  /  (route ( row_index+l , 1 )  - 

route ( row_index, 1 ) )  *  ( speed ( row_index  +1)  - 
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speed ( row_index) )  +  speed (row_index) ;  %  interpolate 

speed 
end 
end 


%%  Angle 

newroutetime ( 1 , 8 )  =  pi  /  2; 

newroutetime ( 1 ,  9 )  =  newroutetime ( 1 , 7 )  /pi  *  180; 

for  i=2 : length (newroutetime ( :  ,  1 ) ) -1 

newroutetime ( i , 5 )  =  newroutetime ( i+1 ,  2 )  - 
newroutetime ( i-1 , 2 )  ;  %  delta  x 

newroutetime (i, 6)  =  newroutetime ( i+1 , 3 )  - 
newroutetime (i-1, 3) ;  %  delta  y 
newroutetime ( i , 7 )  = 

atan2 (newroutetime (i, 6) , newroutetime (i, 5 )) ;  %  heading 
with  horizontal  as  0  rad 

if  newroutetime (i, 7 )  >  pi/2  &&  newroutetime ( i , 7 )  <  pi 

%  heading  with  vertical  as  0  rad 

newroutetime ( i , 8 )  =  2  *  pi  +  (  pi  /  2  - 

newroutetime (i, 7) ) ; 
else 

newroutetime ( i , 8 )  =  pi  /  2  -  newroutetime ( i , 7 ) ; 
end 

newroutetime (i, 9)  =  newroutetime ( i , 8 )  /pi  *  180; 
end 

%%  Rudder 

for  i=l:length (newroutetime ( : , 1 ) ) 

row_index  =  find (newroutetime ( i, 2 ) <route (:, 1 ), 1 ) -1 ; 
if  row_index  <  length ( route (:, 1 ) ) 
newroutetime ( i , 1 0 )  =  (newroutetime (i, 2 )  - 

route (row_index, 1) )  /  (route ( row_index+l , 1 )  - 

route (row_index, 1) )  *  (rudder (row_index  +1)  - 

rudder (row_index) )  +  rudder (row_index) ;  %  interpolate 

speed 
end 
end 

for  i=l : length (newroute (:, 1 ) ) 

row_index  =  find (newroute ( i , 1 ) <route (:, 1 ), 1 ) -1 ; 
if  row_index  <  length ( route (:, 1 ) ) 
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newroute (i, 5)  =  (newroute ( i , 1 )  -  route ( row_index, 1 ) ) 

/  (route (row_index+l , 1 )  -  route ( row_index, 1 ) )  * 

( rudder ( row_index  +  1)  -  rudder ( row_index) )  + 

rudder (row_index) ;  %  interpolate  speed 
end 
end 

TSPI  =  [newroutetime (:, 1 )  newroutet ime ( : , 2 ) 
newroutetime ( : , 3)  newroutetime ( : , 4) . * 
sin (newroutetime ( : , 8 ) )  newroutetime ( : , 4 ) . * 
cos (newroutetime ( : , 8) ) ] ; 
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APPENDIX  F.  MATLAB  CODE  FOR  ROUTE  ANALYSIS 


clear  all 
clc 

%%  Generate  X  Y  Z  from  GPX 

route  =  gpxread ( ' Carmel_Valley . gpx ' ,  ' FeatureType ' , 

' track ' ) ; 

%%  Defining  WGS-84  Geodetic  Constants 

a  =  6378137.0;  %  WGS-84  Earth  semimajor  axis  (m) 

b  =  6356752.3142;  %  WGS-84  Earth  semiminor  axis  (m) 

f  =  (a-b) /a;  %  Ellipsoid  Flatness 

e_sq  =  f  *  (2  -  f ) ;  %  Square  of  Eccentricity 

epsilon  =  e_sq  /  (1  -  e_sq) ; 

h  =  0 ; 

routelLat=route . Latitude ' ;  %  extracting  latitude,  phi 
routelLon=route . Longitude ' ;  %  extracting  longitude, 
lambda 

routelLat2rad=route . Latitude (pi/180) ;  %  extracting 
latitude,  phi 

routelLon2rad=route . Longitude (pi/180 ) ;  %  extracting 
longitude,  lambda 

routel= [ routelLat 2rad  routelLon2rad] ; 
for  i=l : length (routel ) 

routel(i,3)  =  sin (routel (i, 1) ) ;  %  sin  (lat) 
routel  (i, 4)  =  a  /  sqrt(l-e_sq  *  routel  (i, 3)  * 

routel (i, 3) ) ;  %  N,  distance  from  the  surface  to  the  z- 
axis  along  the  ellipsoid  normal 

routel (i, 5)  =  (h  +  routel (i, 4))  *  cos (routel (i, 1 ) )  * 

cos (routel (i, 2 )) ;  %  converting  to  ECEF  x-coordinates 

routel  (i, 6)  =  (h  +  routel (i, 4))  *  cos ( routel ( i , 1 ) )  * 

sin (routel (i, 2 )) ;  %  converting  to  ECEF  y-coordinates 

routel (i, 7)  =  (h  +  (1  -  e_sq)  *  routel  (i, 4))  * 

sin (routel (i, 1 )) ;  %  converting  to  ECEF  z-coordinates 

routel (i, 8)  =  routel (i, 5)  -  routel  (1,5);  %  delta  x 
routel (i, 9)  =  routel (i, 6)  -  routel (1,6);  %  delta  y 
routel (i, 10)  =  routel (i, 7)  -  routel (1,7);  %  delta  z 

routel(i,ll)  =  -sin  ( routel  ( 1 , 2 ) )  *  routel (i, 8)  + 

cos  (routel ( 1 , 2 ) )  *  routel (i, 9);  %  Easting 
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routel(i,12)  =  -cos  ( routel  ( 1 , 2 ) )  *  sin (routel (1, 1) )  * 

routel(i,8)  -  sin (routel ( 1 , 1 ) )  *  sin (routel ( 1 , 2 ) )  * 

routel  (i, 9)  +  cos  (routel ( 1 , 1 ) )  *  routel (i, 10 ) ;  % 

Northing 

routel (i, 13)  =  cos (routel  (1, 1) )  *  cos (routel  (1, 2) )  * 

routel  (i, 8)  +  cos (routel ( 1 , 1 ) )  *  sin  (routel  ( 1 , 2 ) )  * 
routel  (i, 9)  +  sin  (routel ( 1 , 1 ) )  *  routel ( i, 10 ) ;  %  Up 
end 

%%  Convert  Intervals  of  10m 

for  i=l : find (routel (:, 11 ) ==routel (end, 11 )) -1  %finding 

difference  in  x  and  y  for  points  i  and  i+1 

difflx  (i) =routel  (i  +  1, 11) -routel  (i, 11) ; 

dif f ly ( i ) =routel ( i+1 , 12 ) -routel ( i, 12 ) ; 

end 

k=2;  %initiatise  k=2 

carryover ( 1 ) =0 ;  %  initiatise  carryover=0 
addx=0;  addy=0; 

newroute= [routel (1, 11)  routel (1,12)]; 
for  j  =  1 : i 

distance (j) = (difflx (j) A2+diffly (j) A2 ) A0 . 5;  % computing 
distance  between  points  j  and  j+1 
if  carryover (j)  >  0 

if  (carryover ( j ) +distance ( j ) ) >=10 

for  count =1 : floor ( (carryover (j) +di stance (j) ) /10) 
if  count==l 

newroute (k, 1) =newroute (k- 

1. 1)  +difflx ( j) /distance ( j) * (10-carryover ( j) ) +addx; 

newroute (k, 2) =newroute (k- 

1. 2)  +diffly ( j) /distance ( j) * (10-carryover ( j) ) faddy; 

k=k+l ; 
else 

newroute (k, 1) =newroute (k- 

1.1)  fdifflx ( j) /distance ( j) *10; 

newroute (k, 2) =newroute (k- 

1. 2)  fdiffly (j) /distance (j) *10; 

k=k+l; 

end 

end 

carryover ( j+1 ) =di stance ( j ) - (10-carryover ( j ) ) - 
(count-1) *10; 

addx=dif f lx ( j ) / distance ( j ) * carryover ( j+1 ) ; 
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addy=dif f ly ( j )  / distance ( j ) * carryover ( j+1 ) ; 
else 

carryover ( j+1 ) =carryover ( j ) +di stance ( j ) ; 
addx=addx+dif f lx ( j ) ; 
addy=addy+dif f ly ( j ) ; 
end 
else 

if  distance ( j ) >=1 0  %  check  if  distance  to  next 
point  is  >  10m 

for  count =1 : floor (distance ( j ) /10) 
newroute (k, 1) =newroute (k- 

1.1)  +difflx ( j) /distance ( j) *10; 

newroute (k, 2) =newroute (k- 

1.2)  +diffly ( j) /distance ( j) *10; 

k=k+l; 

end 

carryover (j+1) =di stance ( j ) - (count) *10; 
addx=dif f lx ( j ) / distance ( j ) * carryover ( j+1 ) ; 
addy=dif f ly ( j ) / distance ( j ) * carryover ( j+1 ) ; 
elseif  distance ( j ) <10  %  check  if  distance  to  next 
point  is  <  10m 

carryover ( j+1 ) =di stance ( j ) ; 
addx=dif f lx ( j ) ; 
addy=dif f ly ( j ) ; 
end 
end 
end 


%%  Computing  UP  Coordinates  from  Northing 
rindex=0 ; 

newroutedist ( : , 1 )  =  newroute (:, 1 ) ; 
newroutedist ( : , 2 )  =  newroute (:, 2 ) ; 

for  i=l : length (newroutedist )  %  Computing  UP  Coordinates 

from  Northing  for  10m  interval 

rindex=f ind (newroutedist (i, 2 ) >=routel ( : , 12 ) , 1 ) ; 
newroutedist ( i , 3 ) = (newroutedist ( i, 2 ) - 
routel (rindex, 12) ) / (routel (rindex+1, 12) - 
routel (rindex, 12 ) )  *  ( routel ( rindex+1 , 13 ) - 

routel (rindex, 13 ) )  +  routel (rindex, 13 ) ; 
end 
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%%  Heading  Computation  for  Distance-Based 

for  i=2 : length (newroutedist ) -1 

newroutedist (i, 4 )  =  newroutedist ( i+1 , 1 )  - 

newroutedist (i-1, 1 ) ;  %  delta  x 

newroutedist (i, 5)  =  newroutedist ( i+1 , 2 )  - 
newroutedist ( i-1 , 2 ) ;  %  delta  y 

newroutedist (i, 6)  = 

atan2 (newroutedist (i, 5) , newroutedist (i,  4 ))  ;  %  heading 
with  horizontal  as  0  rad 

if  newroutedist (i, 6)  >  pi/2  &&  newroutedist (i, 6)  <  pi 

%  heading  with  vertical  as  0  rad 

newroutedist (i, 7 )  =  2  *  pi  +  (  pi  /  2  - 

newroutedist (i, 6) ) ; 
else 

newroutedist (i, 7 )  =  pi  /  2  -  newroutedist (i, 6) ; 
end 
end 

for  i=2 : length (newroutedist ) -2 

if  newroutedist ( i+1 , 7 )  -  newroutedist ( i, 7 )  >=  pi  && 

newroutedist ( i+1 , 7 )  -  newroutedist ( i , 7 )  <=  3  *  pi  % 

change  in  heading 

newroutedist (i, 8 )  =  newroutedist ( i+1 , 7 )  - 

newroutedist (i, 7 )  -  2  *  pi; 

elseif  newroutedist ( i+1 ,  7 )  -  newroutedist ( i , 7 )  <=  -pi 

&&  newroutedist ( i+1 ,  7 )  -  newroutedist ( i , 7 )  >=  -3  *  pi 

newroutedist (i, 8 )  =  newroutedist ( i+1 , 7 )  - 

newroutedist ( i , 7 )  +  2  *  pi; 

else 

newroutedist (i, 8 )  =  newroutedist ( i+1 , 7 )  - 

newroutedist (i, 7) ; 

end 

end 

if  std (newroutedist (2 : length (newroutedist) -2 ,  8) ) /pi*180 
>=  0  &&  std (newroutedist (2 : length (newroutedist ) - 
2,8))/pi*180  <=  1  %  classify  road  type 
roadtype  =  1; 

elseif  std (newroutedist (2 : length (newroutedist) - 
2,8)) /pi*180  >  1  && 

std (newroutedist ( 2 : length (newroutedist )-2,8))/pi*180  <= 
4 . 5 
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roadtype  =  2; 

elseif  std (newroutedist (2 : length (newroutedist ) - 

2, 8) ) /pi*180  >  4.5  && 

std (newroutedist (2 : length (newroutedist) -2,8))/pi*180  < 

9.5 

roadtype  =  3; 

elseif  std (newroutedist (2 : length (newroutedist) - 

2, 8) ) /pi*180  >  9.5  && 

std (newroutedist (2 : length (newroutedist) -2,8))/pi*180  < 

14.5 

roadtype  =  4; 

elseif  std (newroutedist (2 : length (newroutedist) - 

2, 8) ) /pi*180  >  14.5  && 

std (newroutedist (2 : length (newroutedist) -2,8))/pi*180  < 

19.5 

roadtype  =  5; 

elseif  std (newroutedist (2 : length (newroutedist) - 

2, 8) ) /pi*180  >  19.5  && 

std (newroutedist (2 : length (newroutedist) -2,8))/pi*180  > 

30.5 

roadtype  =  6; 
else 

roadtype  =  7; 
end 

fprintf ( ' \nThe  sigma  change  in  heading  is  %g  degrees, 
therefore  it  is  road  type  %d.\n', 

std (newroutedist (2 : length (newroutedist) -2,  8) ) /pi*180, 
roadtype) 
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APPENDIX  G.  VBA  CODE  FOR  ROUTE  ANALYSIS 


Sub  Main() 

With  ThisWorkbook 
.Sheets("Workspace").  Visible  =  True 
Call  ImportGPX 
Call  CopyColumnToWorkbook 
Call  FillDownXY 
Call  convertm 
Call  HeadingXY 
Call  Chart 
Call  RemoveGPX 

'.Sheets("Workspace"). Visible  =  xlSheetVeryHidden 
End  With 
End  Sub 

Sub  lmportGPX() '  Import  GPX  data  into  worksheet 
With  ThisWorkbook 

. Sheets. Add(After:=.Sheets(. Sheets. Count)). Name  =  "GPX  data" 

End  With 

Application. DisplayAlerts  =  False 

ActiveWorkbook.Xmllmport  Url:="C:\Weapon\GPXfile\Route.gpx", 

lmportMap:=Nothing,  _ 

Overwrite:=True,  Destination:=Worksheets("GPX  data").Range("$A$1 ") 
Application. DisplayAlerts  =  True 

End  Sub 

Sub  CopyColumnToWorkbook() 

Dim  rl  As  Range 
With  ThisWorkbook 
.SheetsfGPX  data"). Select 
Set  rl  =  Range("L:M") 

.Sheets("Workspace"). Select 

.Sheets("Workspace").Range("A4:Z1 000000"). ClearContents 
rl  .Copy  RangefAl ") 

End  With 
End  Sub 


Sub  FillDownXY() 

Dim  strFormulas1(1  To  2),  strFormulas2(1  To  3)  As  Variant 
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Dim  LastRow  As  Long 

Dim  a,  b,  f,  e_sq,  epsilon,  h  As  Long 

a  =  6378137  'WGS-84  Earth  semimajor  axis  (m) 

b  =  6356752.3142  'WGS-84  Earth  semiminor  axis  (m) 

f  =  (a  -  b)  /  a  'Ellipsoid  Flatness 
e_sq  =  f  *  (2  -  f)  'Square  of  Eccentricity 
epsilon  =  e_sq  /  (1  -  e_sq) 

h  =  0 

LastRow  =  RangefA"  &  Rows.Count).End(xlUp).Row 

With  ThisWorkbook.SheetsfWorkspace") 
strFormulasl(l)  =  "=A4/180*pi()" 
strFormulas1(2)  =  "=B4/180*pi()" 

.Range("C4:D4"). Formula  =  strFormulasI 
.Range("C4:D"  &  LastRow).  Fill  Down 

For  j  =  1  To  LastRow  -  3 

.Range("E"  &  j  +  3). Value  =  a  /  ((1  -  e_sq  *  Sin(.Range("C"  &  j  +  3). Value)  * 
Sin(.Range("C"  &  j  +  3). Value))  A  0.5)  'N,  distance  from  the  surface  to  the  z-axis 
along  the  ellipsoid  normal 

.RangefF"  &  j  +  3). Value  =  (h  +  .RangefE"  &  j  +  3). Value)  *  Cos(.Range("C" 
&  j  +  3). Value)  *  Cos(.Range("D"  &  j  +  3). Value)  'converting  to  ECEF  x- 
coordi  nates 

.Range("G"  &  j  +  3). Value  =  (h  +  .RangefE"  &  j  +  3). Value)  *  Cos(.Range("C" 
&  j  +  3). Value)  *  Sin(.Range("D"  &  j  +  3). Value)  'converting  to  ECEF  y- 
coordi  nates 

.RangefH"  &  j  +  3). Value  =  (h  +  (1  -  e_sq)  *  .RangefE"  &  j  +  3). Value)  * 
Sin(.Range("C"  &  j  +  3). Value)  'converting  to  ECEF  z-coordinates 

.Rangefl"  &  j  +  3).Value  =  .RangefF"  &  j  +  3).Value  -  .RangefF"  &  1  + 
3).Value  'delta  x 

.RangefJ"  &  j  +  3).Value  =  .RangefG"  &  j  +  3).Value  -  .RangefG"  &  1  + 
3).Value  'delta  y 

.RangefK"  &  j  +  3).Value  =  .RangefH"  &  j  +  3).Value  -  .RangefH"  &  1  + 
3).Value  'delta  z 

.RangefL"  &  j  +  3).Value  =  -Sin(.Range("D"  &  1  +  3).Value)  *  .Rangefl"  &  j  + 
3).Value  +  Cos(.Range("D"  &  1  +  3). Value)  *  .RangefJ"  &  j  +  3). Value  'Easting 

.RangefM"  &  j  +  3).Value  =  -Cos(.Range("D"  &  1  +  3). Value)  * 

Sin(.Range("C"  &  1  +  3).Value)  *  .Rangefl"  &  j  +  3).Value  -  Sin(.Range("C"  &  1  + 
3).Value)  *  Sin(.Range("D"  &  1  +  3). Value)  *  .RangefJ"  &  j  +  3). Value  + 
Cos(.Range("C"  &  1  +  3).Value)  *  .RangefK"  &  j  +  3).Value  'Northing 
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.RangefN"  &  j  +  3). Value  =  Cos(.Range("C"  &  1  +  3). Value)  * 

Cos(.Range("D"  &  1  +  3). Value)  *  .Rangefl"  &  j  +  3). Value  +  Cos(.Range("C"  &  1 
+  3). Value)  *  Sin(.Range("D"  &  1  +  3). Value)  *  .RangefJ"  &  j  +  3). Value  + 
Sin(.Range("C"  &  1  +  3).Value)  *  .Range("K"  &  j  +  3).Value  'Up 
Next  j 

strFormulas2(1)  =  "=L5-L4" 
strFormulas2(2)  =  "=M5-M4" 
strFormulas2(3)  =  "=SQRT(04A2+P4A2)M 

.Range("04:Q4"). Formula  =  strFormulas2 
.Range("04:Q"  &  LastRow  -  1).FillDown 
End  With 

End  Sub 


Sub  convertm() 

Dim  i,  j,  k,  Count  As  Integer 

Dim  carryover,  addx,  addy,  distance  As  Long 

addx  =  0 
addy  =  0 
k  =  2 

i  =  Range("P"  &  Rows. Count). End(xlUp). Row  -  3 

With  ThisWorkbook.SheetsfWorkspace") 

.Range("R4"). Value  =  0 
.Range("S4"). Value  =  .RangefL"  &  4). Value 
.Range("T4"). Value  =  .RangefM"  &  4). Value 
For j  =  1  To  i 

If  .RangefR"  &  j  +  3). Value  >  0  Then 
If  .RangefR"  &  j  +  3). Value  +  .RangefQ"  &  j  +  3). Value  >=  10  Then 
For  Count  =  1  To  Application.WorksheetFunction.Floor((. RangefR"  &  j  + 
3).Value  +  .RangefQ"  &  j  +  3). Value)  / 10,  1) 

If  Count  =  1  Then 

.RangefS"  &  k  +  3). Value  =  .RangefS"  &  k  +  3  -  1). Value  +  .RangefO" 
&  j  +  3).Value  /  .RangefQ"  &  j  +  3). Value  *  (10  -  .RangefR"  &  j  +  3). Value)  + 
addx 

.RangefT"  &  k  +  3).Value  =  .RangefT"  &  k  +  3  -  1).Value  +  .RangefP" 
&  j  +  3).Value  /  .RangefQ"  &  j  +  3). Value  *  (10  -  .RangefR"  &  j  +  3). Value)  + 
addy 

k  =  k  +  1 
Else 
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.RangefS"  &  k  +  3). Value  =  .Range("S"  &  k  +  3  -  1).Value  +  .Range("0" 
&  j  +  3). Value  /  .RangefQ"  &  j  +  3).Value  *  10 

.RangefT"  &  k  +  3).Value  =  .Range("T"  &  k  +  3  -  1).Value  +  .Range("P" 
&  j  +  3). Value  /  .RangefQ"  &  j  +  3). Value  *  10 
k  =  k  +  1 
End  If 
Next  Count 

.RangefR"  &  j  +  3  +  1). Value  =  .RangefQ"  &  j  +  3). Value  -  (10  - 
.RangefR"  &  j  +  3).Value)  -  (Count  -  2)  *  1 0 

addx  =  .RangefO"  &  j  +  3).Value  /  .Range("Q"  &  j  +  3).Value  *  .RangefR" 
&  j  +  3  +  1).Value 

addy  =  .RangefP"  &  j  +  3).Value  /  .RangefQ"  &  j  +  3).Value  *  .RangefR" 
&  j  +  3  +  1).Value 
Else 

.RangefR"  &  j  +  3  +  1  ).Value  =  .RangefQ"  &  j  +  3).Value  +  .RangefR"  &  j 
+  3). Value 

addx  =  .RangefO"  &  j  +  3). Value  +  addx 
addy  =  .RangefP"  &  j  +  3). Value  +  addy 
End  If 
Else 

If  .RangefQ"  &  j  +  3). Value  >=  10  Then 
For  Count  =  1  To  Application.WorksheetFunction.Floorf. RangefQ"  &  j  + 
3).  Value)/ 10,  1) 

.RangefS"  &  k  +  3).Value  =  .RangefS"  &  k  +  3  -  1).Value  +  .RangefO" 
&  j  +  3). Value  /  .RangefQ"  &  j  +  3). Value  *  10 

.RangefT"  &  k  +  3).Value  =  .RangefT"  &  k  +  3  -  1).Value  +  .RangefP"  & 
j  +  3). Value  /  .RangefQ"  &  j  +  3).Value  *  10 
k  =  k  +  1 
Next  Count 

.RangefR"  &  j  +  3  +  1).Value  =  .RangefQ"  &  j  +  3).Value  -  (Count  - 1)  *  10 
addx  =  .RangefO"  &  j  +  3). Value  /  .RangefQ"  &  j  +  3). Value  *  .RangefR" 
&  j  +  3  +  1).Value 

addy  =  .RangefP"  &  j  +  3).Value  /  .RangefQ"  &  j  +  3).Value  *  .RangefR" 
&  j  +  3  +  1).Value 
Else 

If  .RangefQ"  &  j  +  3). Value  <  10  Then 
.RangefR"  &  j  +  3  +  1 ). Value  =  .RangefQ"  &  j  +  3). Value 
addx  =  .RangefO"  &  j  +  3).Value 
addy  =  .RangefP"  &  j  +  3).Value 
End  If 
End  If 
End  If 
Next  j 
End  With 
End  Sub 
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Sub  Chart() 

Dim  myChtObj  As  ChartObject 
Dim  rngChtData  As  Range 
Dim  rngChtXVal  As  Range 
Dim  rngRefData  As  Range 
Dim  rngRefXVal  As  Range 
Dim  RngToCover  As  Range 
Dim  LastRow  As  Long 

With  ThisWorkbook 
SheetsfWorkspace").  Activate 
LastRow  =  RangefT"  &  Rows. Count). End(xlUp). Row 

Set  rngChtData  =  ActiveSheet.Range("T4:T"  &  LastRow) 
Set  rngChtXVal  =  ActiveSheet.Range("S4:S"  &  LastRow) 
Set  rngRefData  =  ActiveSheet.Range("AA4:AA5") 

Set  rngRefXVal  =  ActiveSheet.Range("Z4:Z5") 
Sheets("Main").  Activate 

Set  myChtObj  =  ActiveSheet.ChartObjects.Add  _ 
(Left:=250,  Width:=375,  Top:=75,  Height:=225) 

With  myChtObj. Chart 

.ChartType  =  xIXYScatterLines 

'  remove  extra  series 
Do  Until  .SeriesCollection. Count  =  0 
.SeriesCollection. Delete 
Loop 

'  add  series  from  selected  range,  column  by  column 
With  . SeriesCollection. NewSeries 
.Values  =  rngChtData 
.XValues  =  rngChtXVal 
.Name  =  "Route  Plot" 

End  With 

With  .SeriesCollection. NewSeries 
.Values  =  rngRefData 
.XValues  =  rngRefXVal 
.Name  =  "1km  Reference" 

End  With 

Set  RngToCover  =  ActiveSheet.Range("G8:023") 
myChtObj. Height  =  RngToCover.Height '  resize 
myChtObj. Width  =  RngToCover. Width  '  resize 
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myChtObj.Top  =  RngToCover.Top  '  reposition 
myChtObj.Left  =  RngToCover.Left  '  reposition 
End  With 
End  With 
End  Sub 

Sub  HeadingXY() 

Dim  I  As  Integer 

Dim  strFormulas1(1  To  3),  strFormulas2(1  To  1),  strFormulas3(1  To  1), 
strFormulas4(1  To  1),  strFormulas5(1  To  1),  strFormulas6(1  To  1)  As  Variant 
Dim  LastRow  As  Long 

With  ThisWorkbook 
.SheetsfWorkspace").  Select 
End  With 

LastRow  =  Range("T"  &  Rows. Count). End(xlUp). Row 

With  ThisWorkbook.SheetsfWorkspace") 
strFormulasl(l)  =  "=S6-S4" 
strFormulas1(2)  =  "=T6-T4" 
strFormulas1(3)  =  "=atan2(U5,V5)/pi()*180" 

.Range("U5:W5"). Formula  =  strFormulasI 
.Range("U5:W"  &  LastRow  -  1).FillDown 
For  I  =  1  To  LastRow  -  5 

If  .Range("W"  &  I  +  4). Value  >  90  And  .Range("W"  &  I  +  4). Value  <  180  Then 
.RangefX"  &  I  +  4).Value  =  360  +  (90  -  .RangefW"  &  I  +  4).Value) 

Else 

.RangefX"  &  I  +  4).Value  =  90  -  .RangefW"  &  I  +  4).Value 
End  If 
Next  I 

strFormulas2(1)  =  "=IF(AND(X6-X5<=540,X6-X5>=180),X6-X5- 

360,IF(AND(X6-X5>=-540,X6-X5<=-1 80),X6-X5+360,  X6-X5))" 

.Range("Y5"). Formula  =  strFormulas2 
.Range("Y5:Y"  &  LastRow  -  2).FillDown 
.RangefV"  &  LastRow  +  2). Value  =  "Mean  Heading" 

.RangefW"  &  LastRow  +  2).Value  =  "Sigma  Heading" 

.RangefX"  &  LastRow  +  2). Value  =  "Mean  Change  Heading" 

.RangefY"  &  LastRow  +  2). Value  =  "Sigma  Change  Heading" 

.RangefY"  &  LastRow  +  4). Value  =  "Road  Type" 
strFormulas3(1)  =  "=AVERAGE(X4:X"  &  LastRow  - 1  &  ")" 
strFormulas4(1)  =  "=STDEVA(X4:X"  &  LastRow  - 1  &  ")" 
strFormulas5(1)  =  "=AVERAGE(Y5:Y"  &  LastRow  - 1  &  ")" 
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strFormulas6(1)  =  "=STDEVA(Y5:Y"  &  LastRow  - 1  &  ")" 

.RangefV"  &  LastRow  +  3). Formula  =  strFormulas3 
.RangefW"  &  LastRow  +  3). Formula  =  strFormulas4 
.RangefX"  &  LastRow  +  3). Formula  =  strFormulas5 
.RangefY"  &  LastRow  +  3). Formula  =  strFormulas6 

If  .Range("Y"  &  LastRow  +  3). Value  >=  0  And  .RangefY"  &  LastRow  + 
3).Value  <=  1  Then 

.RangefY"  &  LastRow  +  5). Value  =  1 

Elself  .RangefY"  &  LastRow  +  3). Value  >  1  And  .RangefY"  &  LastRow  + 
3).Value  <=  4.5  Then 

.RangefY"  &  LastRow  +  5). Value  =  2 

Elself  .RangefY"  &  LastRow  +  3). Value  >  4.5  And  .RangefY"  &  LastRow  + 
3).Value  <=  9.5  Then 

.RangefY"  &  LastRow  +  5). Value  =  3 

Elself  .RangefY"  &  LastRow  +  3). Value  >  9.5  And  .RangefY"  &  LastRow  + 
3).Value  <=  14.5  Then 

.RangefY"  &  LastRow  +  5). Value  =  4 

Elself  .RangefY"  &  LastRow  +  3). Value  >  14.5  And  .RangefY"  &  LastRow  + 
3). Value  <=  19.5  Then 

.RangefY"  &  LastRow  +  5). Value  =  5 

Elself  .RangefY"  &  LastRow  +  3). Value  >  19.5  And  .RangefY"  &  LastRow  + 
3). Value  >=  30.5  Then 

.RangefY"  &  LastRow  +  5). Value  =  6 
Else 

.RangefY"  &  LastRow  +  5). Value  =  "ERROR!!" 

End  If 

If  .RangefV"  &  LastRow  +  3). Value  >=  0  And  .RangefV"  &  LastRow  + 
3). Value  <=  90  Then 

.Range("Z4"). Value  =  .Range("S4"). Value 
.Range("AA4"). Value  =  .Range("T4").Value  -  2000 
.Range("Z5").  Value  =  .Range("S4").  Value  +  1000 
.Range("AA5"). Value  =  .Range("T4").Value  -  2000 

Elself  .RangefV"  &  LastRow  +  3).Value  >  90  And  .RangefV"  &  LastRow  + 
3).Value  <=  180  Then 

.Range("Z4"). Value  =  .Range("S4").Value 
.Range("AA4").  Value  =  .Range("T4").Value  +  2000 
.Range("Z5").  Value  =  .Range("S4").  Value  +  1000 
.Range("AA5").  Value  =  .Range("T4").Value  +  2000 

Elself  .RangefV"  &  LastRow  +  3). Value  <  0  And  .RangefV"  &  LastRow  + 
3).Value  >=  -90  Then 

.Range("Z4"). Value  =  .Range("S4"). Value  -  1000 
.Range("AA4"). Value  =  .Range("T4").Value  -  2000 
.Range("Z5"). Value  =  .Range("S4").Value 
.Range("AA5"). Value  =  .Range("T4").Value  -  2000 
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Else 

.Range("Z4"). Value  =  .Range("S4"). Value  -  1000 
.Range("AA4").  Value  =  .Range("T4").Value  +  2000 
.Range("Z5"). Value  =  .Range("S4").Value 
.Range("AA5").  Value  =  .Range("T4").Value  +  2000 
End  If 
End  With 

Worksheets("Workspace").Range("Y"  &  LastRow  + 
Worksheets("Main").Range("F5") 
Worksheets("Main").Range("F5").lnterior.Colorlndex  =  6 
Worksheets("Main").Range("F5").  Font.  Name  =  "Calibri" 
Worksheets("Main").Range("F5").Font.Size  =  20 
Worksheets("Main").Range("F5").Font.Bold  =  True 
Worksheets("Main").Range("F5").HorizontalAlignment  =  xlCenter 
End  Sub 


Sub  RemoveGPX() 

With  ThisWorkbook 

Application. DisplayAlerts  =  False 
.SheetsfGPX  data"). Delete 
Application. DisplayAlerts  =  True 
\Sheets("Workspace"). Visible  =  xlSheetVeryHidden 
.Sheets("Main").  Select 

End  With 
End  Sub 


5). Copy 
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